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I .  INTRODUCTION 


The  finite  element  method  has  become  a  popular  numerical 
tool  in  the  analysis  of  fluid  flow  problems  (1,  2,  3,  4,  5). 
Particularly  in  the  regime  of  incompressible  flow  this  method 
has  become  very  competitive  with  the  more  established  finite 
difference  methods  due  to  its  simplicity  of  implementation 
and  generality  of  handling  mixed  boundary  conditions  for  com¬ 
plex  geometries  which  favor  nonuniform  meshes  at  points  of 
singularities.  Accordingly,  the  finite  element  discretization 
process  is  used  herein  to  characterize  the  flow  of  a  polymeric 
melt  under  various  geometric  conditions.  The  particular  approach 
is  the  Galerkin  weighted  residual  formation  of  the  non-symmetric 
integral  equations  (6,2),  with  the  penalty  method  used  for  the 
pressure  term  via  sin  approximation  of  the  incompressible  con¬ 
tinuity  equation  (1,3,7) .  Steady,  two  dimensional  flow  is  treated 
and  viscoelastic  fluid  effects  are  modeled  by  employing  an 
Oldroyd  codeformational  stress  derivative  in  a  modified  Maxwell 
constitutive  equation  (8) . 

The  motivation  for  this  analysis  stems  from  the  development 
of  low  cost,  medium  performance,  plastic  gyroscopic  instruments 
at  the  Charles  Stark  Draper  Laboratory.  With  the  exception  of 
the  momentum  wheel  and  electromagnetic  parts,  a  complete  single 
degree  of  freedom  integrating  gyroscope  has  been  designed  using 
glass  filled  polyphenylene  sulfide  parts.  Performance  goals 
are  in  the  range  of  1  -  10  degrees/hour  drift  rate.  Cost  ad¬ 
vantages  derive  from  elimination  of  precision  secondary  machining 
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of  metallic  components  as  well  as  basic  material  costs.  However, 
the  need  for  uniform  physical  and  mechanical  properties  (e.g., 
dimensional  stability,  thermal  conductivity,  mechanical  compli¬ 
ance)  of  the  parts  to  provide  the  required  performance  after 
possible  long-term  storage  in  unair-conditioned  warehouses 
necessitates  the  correlation  of  the  final  material  state  to 
processing  parameters.  Such  knowledge  will  permit  the  rationale 
selection  of  extrusion  parameters,  post  fabrication  treatments, 
and  subsequent  analysis  of  storage  and  service  environmental 
effects  on  instrument  performance.  Figure  1  shows  a  picture 
of  the  typical  plastic  gyroscope  under  consideration.  Figure  2 
depicts  a  typical  injection  molding  process  for  these  gyroscopic 
parts. 

Roy lance  (9)  has  pointed  out  that  the  information  the 
engineer  is  seeking  in  a  flow  analysis  is  the  location  of  regions 
of  elevated  shear  deformation,  which  can  lead  to  mechanical  de¬ 
gradation  and  higher  residual  stresses,  regions  of  stagnation 
and  recirculation,  at  which  overlong  material  residence  and 
thermal  degradation  might  occur,  and  power  requirements  for  the 
fabrication  process  itself.  Also  of  interest  for  the  gyroscope 
application  are  the  effects  of  the  flow  field  on  the  distribution 
of  the  filler  fibers  which  are  carried  along  by  the  drag  of  the 
fluid.  It  is  possible  that  the  zones  of  filler  depletion  or 
enhancement  which  are  observed  in  molded  parts,  can  be  predicted 
and  controlled  by  evaluation  of  the  calculated  velocity  field. 
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In  the  above  regard  the  numerical  analysis  of  polymer 
melts  can  be  broken  down  into  two  general  categories.  First  is 
the  evaluation  of  the  accuracy  of  the  solution  themselves. 
Calculations  are  made  and  compared  to  known  exact  or  approximate 
analytic  solutions.  Typical  of  these  are  the  pipe/channel  flow, 
drag  flow,  and  die  entry  flow.  By  far  most  of  the  numerical 
studies  have  been  in  this  category.  In  the  second  category  is 
the  application  of  numerical  solutions  to  real  problems.  Only 
three  studies  are  known  to  this  author  which  have  aimed  at  apply¬ 
ing  numerical  results  to  actual  polymer  processing.  The  first 
is  the  work  of  Bigg  (10)  who  used  the  Marker  and  Cell  Finite 
Difference  Scheme  to  specify  preferred  operations  for  the  mixing 
of  polymer  solutions  in  a  single  screw  extruder.  The  second  is 
the  National  Science  Foundation/Industry  supported  work  at  Cornell 
University  (10),  also  using  finite  difference  methods  to  evaluate 
mold  filling  and  control  the  location  and  orientation  of  weld 
lines.  Thirdly  is  the  work  of  Caswell  and  Tanner  (12)  who  effec¬ 
tively  used  the  finite  element  method  to  redesign  wire  coating 
dies  to  eliminate  recirculation. 

The  current  work  falls  into  the  first  category  described 
above,  but  the  intention  of  applying  the  numerical  model,  once 
assessed  for  accuracy  and  utility,  is  kept  firmly  in  mind  and 
discussed  throughout  this  report.  To  conclude  the  introduction, 
it  is  also  necessary  to  describe  how  the  current  analysis  fits 
into  the  completely  general  solution.  In  the  injection  molding 
process,  the  flow  is  non-steady  and  non-isothermal  (but  approxi¬ 
mately  adiabatic  within  the  fluid  boundaries) ,  with  advancing 
free  surfaces  until  the  mold  is  completely  filled.  Upstream  of 


the  flow  front  the  fluid  is  completely  surrounded  by  either 
rigid  boundaries  or  adjacent  fluid.  For  an  incompressible 
fluid,  a  complete  numerical  model  must  therefore  account  for 
unsteady,  non-isothermal,  free  surface  effects.  In  addition, 
the  observance  of  a  finite  recoverable  shear  in  the  rheological 
data  of  polymer  melts  indicates  the  need  to  include  viscoelastic 
effects  in  the  model.  For  unsteady  effects,  since  the  Reynolds 
number  (Re)  of  flow  is  always  much  less  than  unity,  a  good  ap¬ 
proximation  is  achieved  by  ignoring  inertia  and  employing  the 
linear  "creeping"  flow  solution.  The  model  that  we  are  eventually 
striving  for  then  is  an  adiabatic,  viscoelastic  solution  with 
changing  surface  boundaries.  Time  is  included  only  as  temperature 
is  conducted  and  convected  and  es  the  velocity  field  is  perturbed 
by  the  changing  boundary.  The  current  work  investigates  the 
viscoelastic  effects  with  the  simplifying  assumptions  of  two 
dimensional,  steady  state  flow. 

To  this  end,  this  report  contains  a  brief  review  of  the 
finite  element  method,  a  discussion  of  the  viscoelastic  consti¬ 
tutive  models  used  in  the  finite  element  equations,  the  details 
of  the  numerical  schemes  used  in  solving  the  equations,  the 
computer  implementation  of  the  numerical  schemes,  a  discussion 
of  calculations  conducted  for  four  flow  geometries  to  assess 
the  numerical  model,  and  an  evaluation  of  the  application  of 
the  numerical  technique  to  the  gyroscope  fabrication. 


4 


II.  THE  FINITE  ELEMENT  METHOD  IN  FLUID  DYNAMICS 

This  section  is  not  intended  to  be  exhaustive  in  nature, 
but  rather  to  review  some  of  the  more  important  features  of 
the  finite  element  method  employed  in  this  work.  References 
may  be  consulted  for  a  more  thorough  treatment  of  the  methods. 

We  begin  by  repeating  that  the  finite  element  method  is 
an  approximate  method  of  solving  the  differential  equations  of 
boundary  and  initial  value  problems  (1,2).  Field  variables  are 
solved  by  dividing  the  bounded  region  into  subsets  (finite 
elements)  which  themselves  are  governed  by  the  differential 
equations.  By  approximating  the  distribution  of  the  variables 
within  each  finite  element  by  a  trial  function,  the  variables 
at  any  point  in  the  element  can  be  determined  by  a  linear  com¬ 
bination  of  the  variable  at  specified  points  on  the  element 
edges.  These  points  are  called  the  nodes  of  the  element;  the 
variables  at  the  nodes  being  determined  by  solving  linear  alge¬ 
braic  equations  formed  by  assembling  all  of  the  elements  into 
a  matrix  equation  of  order  pqm,  where  p  is  the  number  of  elements, 
q  is  the  number  of  nodes  per  element,  and  m  is  the  number  of 
variables  per  node.  The  coefficients  of  the  variables  in  the 
simultaneous  equations  are  the  integrals  of  the  governing  dif¬ 
ferential  equation  taken  over  the  region  of  the  element  which  is 
bounded  by  that  node. 

Mathematically,  we  write  the  discretization  as: 


5 


with  prescribed  boundary  conditions.  In  equation  II. 1,  F  is 
the  governing  differential  equation,  Q  is  the  entire  region 
and  y ^  is  the  region  of  the  finite  element.  Where  physical 
relations  apply  (such  as  the  virtual  work  principle  in  solid 
mechanics) ,  the  equations  cam  be  formed  in  that  basis.  This 
is  the  approach  used  in  references  1  and  9. 


When  the  differential  equation  is  self-adjoint  (can  be 
written  in  the  form  (py 1 )  '  +  qy  +  f  **  0)  with  appropriate 
boundary  conditions  the  equations  can  be  formed  by  an  abbreviated 
variational  principle  by  merely  multiplying  the  differential 
equation  by  the  variation  of  the  independent  variables,  i.e. 


y"[(py')  ■  +  qy  +  f]6ydYi 


-  «I  - 


(II. 2) 


where  I  is  the  integral  of  the  variational  problem  formed  from 
the  governing  differential  equation.  Of  course,  this  is  merely 
stating  that  the  euler  equation  of  the  variational  principle 
is  identical  to  the  governing  differential  equation  (see  [13]) . 
When  the  equations  are  not  self  adjoint,  or  the  boundary  condi¬ 
tions  are  unsuitable,  an  extremum  principle  can  still  be  found, 
unless  odd  number  derivatives  are  present.  In  that  case,  which 
is  the  situation  with  the  complete  Navier-Stokes  equation  with 
convection,  a  true  extremum  principle  does  not  exist  [14]. 
Formation  of  the  finite  element  equations  by  a  variational 
principle  is  the  Ritz  method.  This  method  is  most  useful  for 
the  "creeping"  flow  solution  of  viscous  fluids  where  the  govern¬ 
ing  differential  equation  is  known  to  be  the  euler  equation  of 
the  proper  extremum  principle  [15]. 
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In  the  case  of  the  complete  Navier-Stokes  equation,  the 
method  of  weighted  residuals  is  used  wherein  the  error  which 
remains  after  substituting  appropriate  trial  functions  into 
the  governing  equation  is  orthogonally  projected  to  a  set  of 
weighting  functions  [2].  By  setting  the  inner  product  of  the 
error  and  the  weighting  function  equal  to  zero,  the  approximate 
differential  equation  is  then  satisfied.  Zienkiewicz  [l]  de¬ 
scribes  the  two  most  popular  methods  of  selecting  the  weighting 
functions  as  the  Galerkin  and  Collocation  methods.  Due  to  its 
generality,  the  Galerkin  method  is  the  most  popular  for  formu¬ 
lating  the  finite  element  equations  for  fluid  flow  problems. 
Selecting  this  method  then  the  element  variables  are  approximated 
by  m 

a  =  I  N.C.  (II. 3) 

j-1  ^  3 

where  a  is  the  field  variable  in  the  element,  Cj  are  the  values 
of  the  variable  at  the  node  points  and  Nj  are  the  set  of  trial 
(basis)  functions  which  satisfy  the  element  boundary  condition. 

When  equation  II. 3  is  substituted  into  the  functional  F  of  equation 
II. 1,  we  obtain  in  general: 

n  f  n  f 

1  I  F (a) dyj  -  Z  I  edy,  t  0  (II. 4) 

i=l  *y.  1  i=l  y.  1 
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where  e  is  the  residual  error  of  the  differential  equation. 
Now  using  the  Galerkin  method  of  forming  the  inner  product 
of  the  error  and  the  trial  functions  we  obtain: 


n  r  m 

If  N.F(I  N.C.)dy,  »  0  (k-l,m)  (11.5) 

i-l<  *  j-1  3  3  1 

Yi 

In  this  manner,  we  form  m  times  n  equations  for  the  determi¬ 
nation  of  the  value  of  variable  a  at  the  points  Cj. 

In  selecting  the  field  variables  to  be  approximated 
F recast  [16]  provides  an  excellent  review  of  the  advantages  and 
disadvantages  of  the  different  formulations.  The  governing 
equations  in  an  eulerian  reference  frame  are  continuity 

V*u^  *  0  (II. 6) 

and  momentum 

p[u,t  +  (u*V)u}“  bQ  -  Vp  +  (II. 7) 

where  in  rectilinear  flow: 

u  is  the  velocity  vector 

b  is  the  body  force  vector 


is  the  deviatoric  stress  vector 


a 


°ii±  *  ai2—  +  “ijS 
°21±  +  °22—  +  °23^ 
0  3i—  +  °32i  +  °33i 


p  is  the  constant  density 
p  is  the  hydrostatic  pressure 

8  S  8 

V  is  the  gradient  operator  ^i  +  j  +  k 

and  the  comma  denotes  differentiation  with  respect  to  time. 

If  the  flow  is  purely  viscous,  the  deviatoric  stresses  can  be 
written  as  explicit  functions  of  the  velocity  gradients  leaving 
only  velocity  and  pressure  as  independent  variables.  If  both 
are  approximated  by  the  Galerkin  method,  the  number  of  unknowns 
is  relatively  high  (i.e.  components  of  velocity  at  each  node 
plus  the  pressure) .  In  addition,  some  of  the  diagonal  terms 
of  the  coefficient  matrix  become  zero  which  limits  the  pivoting 
techniques  generally  used  for  solving  the  equations. 


Two  methods  have  been  devised  for  eliminating  the  pressure. 
For  two  dimensional  flow,  the  stream  function  u  *  i|»,y  and  v  *  -^,x 
is  used  to  satisfy  continuity  and  results  in  the  disappearance 
of  the  pressure  term  when  inserted  into  the  momentum  equation. 
However,  the  application  to  mixed  boundary  value  problems  is 
difficult,  as  shown  by  Tanner  [17].  For  incompressible  problems, 
the  penalty  function  formulation  has  been  developed.  This  method, 
reviewed  in  detail  in  [7],  replaces  the  incompressible  continuity 


equation  by  the  approximation 
p  =  -a  (V*u) 
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(II. 8) 


where  a  is  a  large  positive  number  whose  effect  is  to  "penalize" 
the  error  of  not  satisfying  continuity.  In  reference  [7],  it  is 
shown  that  this  method  converges  to  the  exact  solution  for  "creep 
ing"  flow  and  that  the  selection  of  a  is  determined  from  the 
relation : 

a  =  cy  (II. 9) 

Where  c  is  a  constant  equal  to  107  and  y  is  the  dynamic  viscosity 
Furthermore,  to  avoid  the  trivial  solution  of  u  +  O  as  #  +  « 

(see  equation  II. 10)  the  coefficients  determined  from  evaluation 
of  the  integral  must  be  singular.  This  is  accomplished  by  employ 
ing  reduced  integration  (quadrature  rule  of  lower  order  than 
the  exact  for  a  given  element)  for  the  pressure  term.  The  other 
terms  can  then  be  integrated  at  the  optimum  order  (selective 
reduced  integration)  or  at  the  lower  order  (uniform  reduced  in¬ 
tegration)  .  While  it  is  more  accurate  to  employ  selective 
reduced  integration  (SRI) ,  it  is  usually  more  convenient  to  use 
uniform  reduced  integration  (URI)  in  the  computer  programs. 

Since  it  has  been  shown  that  8  node  quadrilateral  elements 
exhibit  inferior  behavior  to  9  node  elements  even  for  SRI,  it 
is  strongly  recommended  that  when  URI  is  used  the  9  node 
"Lagrange"  isoparametric  element  be  employed  [7]. 

Bercovier  [18]  has  recently  concluded  that  the  reduced 
integration  approach  is  only  valid  for  straight- sided  elements 
(biquadratic)  if  the  governing  equation  is  linear  ("creeping" 
flow)  and  valid  only  for  rectangular  elements  (vice  bilinear 
quadrilaterals)  when  the  equation  is  non-linear  (with  convection) 
Since  most  of  our  work  concerns  linear  systems,  this  is  not 
viewed  as  a  limitation. 
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For  ease  of  implementation#  economy#  and  accuracy#  therefore# 
we  selected  the  penalty  method  with  URI,  9  node  Lagrange  isopara¬ 
metric  elements.  For  comparison,  some  eight  node  "serendipity" 
element  cases  were  run  and  will  be  discussed  in  Section  VI. 

Applying  the  Galerkin  formulation  of  the  finite  element 
equations  we  obtain  the  following  for  two  dimensional#  recti¬ 
linear,  incompressible,  viscous  flow: 

(K  +  K  +  K)  u  +  M-|^u+f  =  0  (II. 10) 


Where 


u  is  a  column  vector  of  the  two  dimensional  velocities 
at  the  node  points, 

N  is  the  matrix  of  trial  (shape)  functions, 

£  =  fy  1T  g  |  dy 


and 


|  *  fr  p  NT  (V-  (s  u)T)Tg  dy 
|  «  /  (mT  g)TamT  g  dy 

y  -  fy  HTP  g  dy 
£  nT  b o  dy  ‘  /  ST  -  d 
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In  the  matrix  definitions  above,  we  used  the  further  identi 


ties: 


D  =  u 

2  0  0 

0  2  0 

,  B  =  L  N  and  m= 

1 

1 

0  0  1 

0 

where 


L  is  the  differential  operator  matrix  for  two  dimensional 


flow 


L 


3x 

o 


o 


3 


1_ 

I  3y  3x  I. 

Also  the  second  term  in  the  expression  for  f  is  the  surface 
traction  on  the  line  element  r  which  results  from  integrating 
the  viscous  stress  term  by  parts.  (Throughout  this  report  a 
single  underline  denotes  a  vector  quantity,  and  a  double  under¬ 
line  denotes  a  matrix  quantity.) 


When  the  inertial  effects  are  comparable  to  the  viscous 
ones,  i.e.,  Reynolds  No.  greater  than  one,  equations  II. 10  are 
non  linear  and  must  be  solved  by  some  iterative  scheme.  A 
discussion  of  these  techniques  will  be  postponed  until  the  non¬ 
linear  viscoelastic  effects  are  added  in  Section  IV. 


Of  course  equation  II. 10  is  the  well  known  "weak"  form 
of  the  Navier-Stokes  governing  differential  equation  which  has 
been  derived  elsewhere  by  the  virtual  work  statement  [l]. 
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III.  VISCOELASTIC  CONSTITUTIVE  MODELS 


The  selection  of  a  viscoelastic  constitutive  model  (the 
rheological  equation  of  state)  for  use  in  the  finite  element 
equations  is  generally  a  compromise  between  the  accuracy  of  the 
model  and  ease  of  implementation.  Because  all  of  the  models  are 
nonlinear  consideration  must  be  given  to  the  relative  effects 
on  the  numerical  convergence  of  the  solutions.  In  this  study, 
two  general  ground  rules  were  used  in  selecting  the  appropriate 
model.  First,  for  the  material  under  consideration  (fiber-filled 
polyphenylene  sulfide) ,  adequate  rheological  or  viscometric  data 
do  not  exist  to  justify  the  use  of  multiple  constant  models,  and 
second  only  a  first  order  effect  on  the  flow  field  was  being 
sought.  Once  success  is  achieved  in  modeling  viscoelasticity, 
rheological  data  can  be  obtained  and  adjustments  to  the  consti¬ 
tutive  model  investigated. 

As  before,  only  essential  elements  for  understanding  the 
behavior  of  the  selected  viscoelastic  model  are  presented  in 
this  report.  For  a  thorough  discussion  of  the  continuum  mechanics 
of  viscoelastic  materials  the  references  can  be  consulted  (19, 

20,  21,  22,  23). 

For  a  fluid  element,  the  resistance  to  deformation  when  a 
force  is  applied  can  be  thought  of  as  a  combination  of  viscous 
and  elastic  stresses.  Modeling  these  as  a  dashpot  and  spring 
respectively  as  shown  in  Figure  3,  we  obtain  the  well-known 
Maxwell  Element  for  fluids.  Using  the  nomenclature  of  Figure  3, 


where  y  is  the  dynamic  viscosity,  G  is  the  shear  modulus  of 
elasticity,  e  is  the  infinitesimal  strain  and  a  is  the  applied 
shear  stress  we  obtain  the  stress-strain  rate  relation: 


Generalizing  to  a  three-dimensional  form,  we  have: 


(III.l) 


£  +  x  It  (£)=  (III. 2) 

where  0  is  the  Cauchy  deviatoric  stress  tensor 
y  is  the  dynamic  viscosity 

X  =  y/G  is  a  time  constant  known  as  the  relaxation  time 
and  d  is  the  rate  of  deformation  tensor  whose  components  are 
defined  as: 


1,2,3) 


(III. 3) 


Equation  III . 2  is  suitable  when  the  rate  of  deformation 
in  the  flow  is  inf initesimually  small.  But  for  general  motion, 
in  which  the  rate  of  deformation  is  not  necessarily  small,  the 
time  derivative  of  the  Cauchy  stress  tensor  violates  two  funda¬ 
mental  requirements  of  any  equation  of  state.  These  requirements 
are  that  the  equation  describes  material  properties  independent 
of  the  frame  of  reference,  and  that  the  behavior  of  any  material 
element  must  depend  only  on  its  previous  deformative  history  and 
not  in  any  way  on  the  state  of  neighboring  elements,  or  on  rigid 
body  translation/rotation.  These  discrepancies  are  corrected  by 
substituting  for  the  time  derivative  of  the  Cauchy  stress  either 
an  Oldroyd  derivative  [8]  (known  as  a  convected  or  a  codeforma- 
tional  derivative)  or  a  Jaumann  derivative  [24]  (known  as  a 


co-rotational  derivative) .  These  modifications  will  be 
discussed  shortly.  Once  the  above  requirements  are  satisfied 
it  only  remains  to  tailor  the  equation  so  as  to  fit  experimental 
observations.  This  is  done  by  introducing  added  parameters 
which  are  multiplied  by  functions  of  the  invariants  of  the  rate 
of  strain  tensor. 

Han  [23]  presents  a  survey  of  the  major  refinements  developed 
for  the  two  invariant  stress  derivatives  along  with  the  material 
properties  they  predict.  A  two  constant  (A,y)  model  using  an 
Oldroyd  derivative  is  known  as  a  White-Metzner  model.  When  the 
Jaumann  derivative  is  used,  the  equation  is  called  a  DeWitt  model. 
As  multiple  parameters  are  added,  the  general  models  are  known 
merely  as  n-order  Oldroyd  models.  Two  other  models  derived  by 
means  somewhat  different  from  the  generalized  Maxwell  element 
are  the  Spriggs  model  which  builds  many  Maxwell  elements  at  the 
molecular  structure  level  and  the  Rivlin  Erickson  fluid  which 
merely  states  that  the  fluid  stress  is  a  function  of  the  invariants 
of  the  gradients  of  displacement,  velocity,  acceleration,  second 
acceleration,  and  so  on. 


Returning  to  the  invariant  stress  derivatives,  we  write  them 
explicitly  for  further  discussion.  For  the  Oldroyd  derivative 
in  contravariant  form  (see  [22]  for  a  discussion  of  covariant 
and  contravariant  tensors)  we  obtain: 


Ki_  8oii 

Tt 


+  u 


3oij  .  !^i  _  „  i 

k  3xk  °kj  ax^  ik  Jx^ 


(III. 4) 


3t  r  “k  kj 

Where  the  range  and  summation  indicial  convention  is  used. 


Similarly  the  Jaumann  derivative  is: 


fZii  =  IZii  +  u 

Pt  3t  uk  9x 


+  wik  aik  +  ujk  °ik 


(III. 5) 


1  /3u i  9ui\ 

where  uk  ^  =  j  l - )  are  the  components  of  the  vorticity 

tensor.  Again,  see  Han  [23]  for  an  excellent  discussion  of  the 
physical  significance  of  the  terms  on  the  right-hand  side  of 
equations  III. 4  and  III. 5. 

We  will  also  have  occasion  to  discuss  further  the  Rivlin- 
Ericksen  fluid  so  we  list  the  general  equation  for  an  incom¬ 
pressible  fluid: 

a  A(1)  +  a2  A^J  +a3  A(2)  +«4  A{2)  +  a5(A(l)A(2)  +  A(2)A(1) 
+  a6  (AU)A(2)  +  A(2)Aa))+a7  (A(2)A(1)  +  A(1)A(2)) 


.2  .  .2  ,  .2  ,2  . 

(1)  A{2)  +  A(2)  +  A(2)  A(l) 


+  0t8  (A 

where  the  a.  are  functions  of  the  invariants  of  A 


(1) 


(III. 6) 
and  Aj2j 


and 


(1) 

Aij  =  2dij 

(2) 

Aij 


(1) 
9Ai  j 
3x. 


(1) 


3uk 


(1) 


3uk 


+  A,'.'  +  A.r  -5 — — 

'k  *3  3xi  lk 

In  passing,  it  is  noted  that  the  preceding  discussion  of 
models  has  focused  on  the  rate  type.  If  equation  III. 2  is  inte¬ 
grated  with  respect  to  time  rheological  equations  of  state  of 
the  integral  type  are  obtained.  While  this  type  proves  useful 
for  some  rheological  investigation,  it  complicates  finite  element 
calculations  by  requiring  a  complete  time  history  of  the  strain 
path  of  all  elements. 
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Finally,  before  we  can  discuss  the  relative  merits  of  the 


models,  we  must  make  some  definitions.  Steady  simple  shear 
flow,  also  known  as  viscometric  flow,  is  defined  by  the  velocity 
field 

u  =  -yy,  v  =  w  =  o  (III. 7) 

where  y  is  a  constant  shear  strain  rate  and 

u  is  the  velocity  normal  to  the  y  axis  of  the  cartesian 
coordinate  system. 

Substituting  equation  III. 7  into  III. 3  we  find  the  rate  of 
deformation  tensor  to  be: 


(III. 8) 


For  viscometric  flow,  viscoelastic  fluids  are  observed  to 
exhibit  three  independent  material  properties,  the  standard 
viscosity,  and  a  first  and  second  normal  stress  function  written 
consecutively  as : 


•  •  .  .  2  •  •  2 

°12  =  "  °22  =  ^i(y)y  '  a22  "  a33  =  ^2(Y>Y 


(III. 9) 


Implicit  in  equations  III. 9  is  the  further  observation  that 
when  a  fluid  behaves  viscoelastically ,  the  material  parameters 
are  not  constant,  but  vary  with  the  rate  of  strain.  This  non- 
newtonian  behavior  is  generally  observed  to  follow  a  power  law 
relation,  written  for  the  viscosity  as: 


U  (Y) 


(III. 10) 


1  +  ( K/y o )  (Y/  2) 

where  y0»  K,  and  n  are  parameters  selected  emperically.  When 
the  exponent  n  is  less  than  zero,  the  viscosity  varies  inversely 
to  the  shear  strain  rate  and  the  fluid  is  termed  shear-thinning, 
when  n  is  greater  than  zero  the  fluid  is  shear  thickening.  Most 
real  fluids  are  shear  thinning.  and  on  the  other  hand  are 
observed  to  increase  exponential lv  with  shear  strain  rate. 

Before  we  continue,  recall  that  equation  III. 10  was  written 
for  simple  shear  flow.  This  equation  is  merely  the  specializa¬ 
tion  of  the  more  commonly  written  general  flow  form: 


y(IId)  = 


1  +  (K/y,,)  <>S  Hd) 


[l-n)/2 


(III. 11) 


where  11^  is  the  second  invariant  of  the  rate  of  strain  tensor 


“a  ’  aijdij' 


(III. 12) 


which  in  two  dimensional  rectilinear  flow  can  be  written  explic¬ 


itly  as  : 


-  «[(£)*+ 


(III. 13) 


We  are  now  prepared  to  make  a  selection  of  the  constitutive 
equation  to  implement  in  the  finite  element  equations.  The 
choices  have  been  narrowed  to  (i)  White-Metzner  (ii)  DeWitt  and 
(iii)  Rivlin-Ericksen  as  generally  representative  of  the  avail¬ 
able  models  (Pipkin  and  Tanner  [25]  present  a  thorough  review 
of  all  the  models  for  viscosmetric  flow) .  Middleman  [26]  has 
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presented  an  excellent  discussion  of  the  correlation  of  the 
properties  predicted  by  the  White-Metzner  and  DeWitt  models 
to  experimental  observations.  In  simple  shear  flow,  the  DeWitt 
model  is  somewhat  superior  because  the  second  normal  stress  func¬ 
tion  is  finite  whereas  the  White-Metzner  model  predicts  that  it 
vanishes.  However,  in  general  flow  fields  the  DeWitt  model 
varies  appreciably  from  reality  while  the  White-Metzner  model 
maintains  consistency.  Since  ip2  is  generally  small,  the  fact 
that  the  White-Metzner  model  predicts  a  zero  value  is  not  con¬ 
sidered  a  major  drawback  by  Middleman  and  we  agree.  Han  [23] 
suggests  that  since  the  Oldroyd  derivative  takes  a  different 
form  if  written  in  terms  of  covariant  or  contravaraint  basis 
vectors  that  it  is  inferior  to  the  Jaumann  derivative.  Since 
the  work  herein  is  conducted  for  a  rectilinear  coordinate 
system,  it  is  felt  that  this  is  less  of  a  penalty  than  the 
cited  deviation  of  the  Jaumann  derivative  model  for  general  flow 
•fields.  Therefore,  the  author  concurs  with  Middleman's  recom¬ 
mendation  that  the  White-Metzner  model  is  preferred  to  the 
DeWitt  model. 

Considering  the  Rivlin-Ericksen  fluid.  Tanner  [17]  notes 
that  for  a  simple  shear  flow  equation  III. 6  reduces  to: 

«lj  *  ^ij1’  *  <*1  +  -  *1  \j2’  <III'11> 

Clearly  equation  III. 6  is  overly  complicated  for  our  initial 
work.  But  sin»a  the  simplification  to  III. 11  presumes  simple 
shear  flow,  it  is  disqualified  as  a  candidate  for  this  effort. 

It  is  interesting  to  note,  however,  that  of  the  three  models 


considered,  the  Rivlin  Ericksen  fluid  alone  permits  the  de- 
viatoric  stress  to  be  written  as  an  explicit  function  of  the 
velocities  and  nth  order  derivatives  of  velocities.  The  ad¬ 
vantages  of  this  fact  will  become  obvious  in  the  next  section 
when  we  discuss  the  formation  and  solution  of  the  complete 
finite  element  eauation. 


Let  us  recapitulate  before  concludina  this  section.  A 
White-Metzner  modified  Maxwell  element  was  selected  for  the 
rheological  ecruation  of  state  because  of  its  ability  to  approxi¬ 
mate  real  viscoelastic  fluid  behavior  while  requiring  onlv  two 
model  parameters.  In  addition,  the  two  parameters  y  and  X  are 
taken  to  be  functions  of  the  second  invariant  of  the  rate  of 
strain  tensor  as  defined  in  equation  III. 11. 


3  3 

For  plane,  steady  flow  where  w  ■  ^  ■  o,  the  nine 

equations  of  III. 2  reduce  to  four  which  are  written  explicitly 


below  with  the  use  of  equation  III. 4. 
3c  3a 


v  +  2 


xx  xx  9u  3u 

3x  v  3y  -  2a  3x  ”  yx  3y 


3u 

57 


(III. 12a) 


a  +  X 
xy 


3a  3a 

xy  ,  „  xy  _  3u 
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3q 
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)• 
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These  equations  are  identical  to  those  used  by  Perera  and  Strauss 
[27]  in  their  finite  difference  formulation  of  similar  problems 
when  account  is  made  of  the  reduction  of  the  four-constant  model 
they  used  vice  the  two  parameter  model  used  herein. 

The  reader  is  reminded  that  the  stresses  in  equations  III. 12 
are  the  deviatoric  ones  and  differ  from  the  complete  stresses  by 
the  hydrostatic  pressure.  Since  the  momentum  equation  always 
expresses  these  two  stresses  separately,  they  are  not  combined 
here  either. 
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IV.  VISCOELASTIC  FINITE  ELEMENT  EQUATIONS 

The  governing  differential  equations  for  an  incompressible 
viscoelastic  fluid  are  as  presented  in  equations  II. 6  and  II. 7. 
Continuity  and  momentum  are  repeated: 

7*u  -  o  (IV. la) 

p  [u,fc  +  (u  •  V)u]  »  ^  -  7p  +  7*o  (IV. lb) 

The  boundary  conditions  of  course  will  be  for  the  independent 
variables  and  gradients  of  these  variables.  However,  for  many 
flow  problems,  it  is  more  convenient  to  specify  the  tractions 
(stresses)  on  some  boundaries  and  the  independent  variables  on 
others.  This  is  the  mixed  boundary  condition  formulation  and 
is  of  course  mandatory  for  finite  element  equations  which  are 
reduced  to  a  set  of  inhomogeneous  linear  algebraic  equations. 
While  specification  of  the  variables  (u,  t|>,  p,  a  depending  on 
the  type  of  equations  used)  at  the  boundaries  is  straight 
forward,  the  specification  of  boundary  tractions  must  be  con¬ 
sistent  with  the  type  of  problem.  For  example ,  Chang  [15] 
discusses  the  difference  in  specifying  the  surface  traction, 
for  a  number  of  flow  cases,  between  a  non-newtonian  viscous 
fluid  and  a  generally  viscoelastic  one.  Understand ing  these 
differences  is  particularly  important  when  a  specific  type  of 
flow  is  prescribed  (e.g.,  fully  developed  entry  flow)  for  an 
assessment  of  the  accuracy  of  the  finite  element  model.  We 
defer  further  comment  on  the  boundary  conditions  until  Section  VI 
when  specific  flow  problems  are  considered. 
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Briefly  reviewing  the  past  work  on  finite  element  modeling 
of  viscoelastic  flow,  it  is  noted  that  no  investigations,  known 
to  the  author,  have  been  conducted  using  the  "penalty"  method 
for  incompressible  fluid  flow.  Tanner  [17 J  and  Caswell  and 
Tanner  [12]  have  used  the  formulation  with  velocities  and  pressure 
as  the  independent  variables,  with  a  Rivlin-Ericksen  fluid  for 
viscometric  flow.  Results  have  been  excellent  for  power  law 
fluids,  but  only  Poiseuille  flow  has  been  considered  for  the 
viscoelastic  case.  Kawahara  and  Takeuchi  [28]  applied  a  mixed 
method  where  the  total  deviatoric  stress  (viscous  and  elastic) 
was  independently  interpolated  along  with  the  velocities  and 
pressure.  The  White-Metzner  constitutive  equation  was  then 
solved  simultaneously  with  the  Navier-Stokes  equation  for  in¬ 
compressible  fluids.  Using  six-node  triangular  elements  in 
plane  flow,  this  gives  rise  to  eighteen  additional  unknowns 
per  element  and  is  felt  to  have  limitations  for  general  problems 
because  of  the  computer  capacity  required  for  large,  complicated 
geometric  problems.  However,  they  did  achieve  good  results  for 
expanding  and  bending  flow  through  channels  for  relaxation  times 
up  to  0.1  seconds. 

In  the  work  most  similar  to  the  current  effort,  Chang  et.  al. 
[15]  solved  the  equations  using  the  White-Metzner  model  with 
velocities  and  pressure  the  field  variables  for  the  finite  ele¬ 
ment  equations.  In  two-dimensional,  steady  state,  convective, 
isothermal  flow,  the  slip  stick  problem  was  solved  for  material 
cases  of  Weissenberg  numbers  up  to  0.2. 


The  Weissenberg  number  is  a  dimensionless  ratio  of  re¬ 
coverable  or  elastic  shear  stress  to  total  shear  stress  in 


steady  flow.  It  is  written 

ws  =  if  (IV. V 

It 

where  X  is  the  relaxation  time  in  seconds, 

U  is  a  characteristic  steady  velocity  in  cm/sec, 
and  L  is  a  characteristic  length  in  cm. 

Han  [23]  presents  rheological  data  for  high  and  low  density 
polyethylene  at  various  shear  strain  rates  (U/L)  at  200°C. 

For  high  density  polyethylene,  the  Ws  varies  from  35  at  0.025 
cm/cm-sec  down  to  0.01  at  100  cm/cm-sec.  On  the  other  hand, 
the  Ws  for  low  density  polyethylene  varies  between  5  at  the  low 
strain  rate  and  0.01  at  the  high  strain  rate.  We  note  that  this 
is  essentially  the  range  of  interest  for  practical  problems 
(0. 0l£WS£35) .  A  major  difficulty  in  the  finite  element  method 
has  been  obtaining  numerical  convergence  for  problems  of  rela¬ 
tively  high  Ws  as  evidenced  in  the  above  review.  It  appears 
that  Chang's  work  has  provided  the  highest  value.  Without  dis¬ 
cussion,  it  is  noted  that  with  this  convergence  problem,  the 
added  numerical  problems  associated  with  evaluation  of  the 
pressure  term  in  the  tangent  stiffness  matrix  for  the  penalty 
method  may  suggest  some  limitations  in  the  future  for  applica¬ 
tion  to  viscoelasticity. 


Now  using  the  Galerkin  formulation  with  the  penalty  method, 
equations  IV. 1  become  for  steady  state 


i  (7* (N  u)T)TN  +  (mT  B)TamT  B)dv 


u- 


adv=o 


(IV. 2) 


Where  all  terms  have  been  defined  in  equation  II. 10,  the 
body  forces  are  assumed  to  be  zero,  and  two-dimensional  recti¬ 
linear  flow  is  treated  so  that  the  plane  stress  vector  a  is: 


a  = 


xx 


yy 


LoxyJ 


(IV. 3) 


We  now  split  the  deviatoric  stess  into  a  viscous  and  elastic 
portion 


v  .  e 
a  =  a  +c 


(IV. 4) 


substitute  into  .equation  IV. 2,  and  apply  Green's  divergence 
theorem  to  obtain 

j  ji ITC  &  +  HT0(5’S*a)T)'rB  +  (mTB)  TdmTa)  dv  j  G  +  yHTLT£edv-^TtdA= 

V  v  (IV.  5) 

where  the  viscous  stress  has  been  written  explicity  as 

/S  A 

a  =  DLNu  =  DBu  (IV. 6) 

and  the  last  term  is  the  traction  on  the  boundary.  Prom  equation 


III. 2  we  can  write 


(av  +  oe)  +  X 


*t~ 


=  2ye 


or  since  a  =  2y£ 

e  &  - 

a  -  -  X-r— r 


(IV. 7) 


(IV. 8) 


where  e  is  the  2D  rate  of  deformation  vector. 
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From  equations  III. 12,  we  see  that  for  steady  state  equation 
IV. 8  is  of  the  following  functional  form 


£®  =  g(£/  £S»  £V*  £@  /  £V#,  u',  x,  y)  (IV. 9) 

Where  the  prime  denotes  differentiation  with  respect  to  x  and  y. 
But  since  £V  is  a  unique  function  of  u'  we  can  further  state 

£e  =  h(£,  u',  u  ,  x,  y,  oe,  ae  ).  (IV. 10) 

Equation  IV. 10  now  makes  equation  IV. 5  not  only  non-linear 
(even  for  creeping  flow) ,  but  inexpressible  in  an  explicit  form. 
The  equation  must,  therefore,  be  solved  simultaneously  with 
equation  IV. 5.  This  is  the  same  point  reached  by  Chang  [15] 
and  Perera  [27].  Let  us  examine  the  method  of  solution  proposed 
in  [15].  Although  convection  was  included  in  that  analysis,  it 
is  easier  to  consider  creeping  flow  (without  loss  of  generality) . 


The  creeping,  viscoelastic  flow  can  be  written  as: 

*  ~  ~  *r  ''tt  a 

S  H  +  5  2  '  H  '  £  '  £  )=  f 


(iv. ll) 


where  the  terms  Ke  are  the  functional  form  of  the  internal  elastic 
forces.  Newton-Raphson  iteration  can  not  be  employed  to  solve 
IV. 11  because  of  the  implicit  dependent  variable  ae.  Instead 
the  common  method  is  to  use  successive  substitution  where  an 
initial  value  of  oe  is  guessed  and  substituted  into  equation 

A 

IV. 10.  Assuming  u  has  first  been  solved  for  the  linear  problem, 

Ke  can  now  be  calculated,  substituted  into  equation  IV. 11  and 

A  A 

a  new  value  of  u  found.  This  new  value  of  u  is  than  used  with 
the  latest  value  of  ae  to  calculate  an  updated  value  of  ce  and 
the  process  is  repeated  until  some  convergence  criterion  is 


es  /  Ag  ^/#s  es— 1  \ 

and  a  *  h(  u  ,  u  ,  u  ,  x,  y,  o  ,  o  I.  (IV. 12) 

The  actual  calculation  on  the  computer  was  performed  at  the  iteration 

es  -s+i 

s+1  by  subtracting  K  from  f  and  solving  K  u  .  Therefore, 

the  computer  equation  is: 

~  es  '"s 

K  A  u  =  f  -  K  -  K  us 

*  Afi+1  ''s 

where  Au  =  u  -  u  .  If  the  convection  non-linearity  is  in¬ 
cluded  the  Picard  substitution  can  be  nested  within  a  Newton- 
Raphson  iteration. 

If  we  momentarily  disregard  the  issue  of  convergence,  the 
only  problem  which  remains  is  the  calculation  of  the  elastic 
stress  gradient  at  the  s-1  iteration.  Chang  [15]  is  completely 
silent  on  this  issue  and  it  is  felt  that  it  was  ignored.  Later 
on,  we  will  discuss  possible  situations  where  this  might  be  valid. 

To  aid  in  the  discussion,  let  us  write  equation  III. 12  in  vector 
form  by  recognizing  a  -  o  .  It  can  be  verified  that  the  equa- 

xy  y x 

tion  becomes: 

a®  =  A[A  £  -  (u  •  V)o]  (IV. 13) 
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Assuming  the  integral  could  be  evaluated  numerically 
could  be  solved  by  the  same  successive  substitution  scheme 
used  for  the  complete  finite  element  equations.  An  initial 
guess  is  made  for  a®x  in  the  integrand  and  the  right-hand  side 
is  solved  for  an  updated  value  of  o®x.  That  value  is  then 
substituted  into  the  integrand  and  the  procedure  repeated  until 
convergence  is  achieved.  Let  us  now  write  IV. 13  in  this  form 


(u*V)a  =  j  (ae  -  A  a)  / 


(IV. 16) 


and  upon  integration  by  taking  the  dot  product  of  both  sides 


with  dA  =  dxi  +  dyj_ 


(u  +  v) 


2*2°  *f  H2"  -  '  “)■' 


(IV. 17) 


While  in  theory,  IV. 17  could  be  solved,  it  is  felt  that  in  a  finite 

element  formulation,  it  would  be  impractical  to  use  such  a  system 

that  requires  an  initial  value  to  be  calculated  at  a  corner  of 

each  element  (ct^  and  separate  integration  of  the  spatial  derivatives, 
i.e. , 

A  x  y  A 

/  x^£e  -  A  0^-dA  =J  X  ^  £e-A  £^dx  +J  x^®-*  £^dxdy 


Due  to  the  difficulties  encountered,  another  method  was 
sought  for  the  solution  of  IV. 13.  If  the  derivative  is  approxi¬ 
mated  by  a  Taylor  series,  then  a  standard  finite  difference 
equation  is  achieved  and  usual  relaxation  methods  can  be  employed 
for  the  solution.  Referring  to  Figure  4  and  using  central  dif¬ 
ferences  we  have  for  the  first  component  of  c® 
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(IV. 18) 

Where 

3o  =  (gi+1^-qi-1^)(yi^+1~vi^"l  )-(oi  '  j+1-oi J'1)  (yi+1  >  j-y1"1 
3x  (  x1+1  '^-x1-1' 3)  (y1 ' ^+1-yi '  j“1  ) -( y1+1 ' 3-y1-1 ' 3)  ( x1  '^-x1 ' 3 

and 

9y  i^j“(xi+1^V”1^j(yx'3+1-y1,3“l)-(yi+l'i-yi"l,j)(x1'3+l-x1'j'l 

(IV. 19) 

Equations  IV. 19  are  derived  in  Appendix  I. 

A  few  words  about  equations  IV. 18  and  IV. 19  are  in  order. 

While  central  differences  are  expected  to  give  higher  order 
accuracy,  Roache  [30]  notes  that  the  numerical  stability  is  much 
poorer  than  backward  (upwind)  differencing  and  for  a  non-uniform 
mesh  (special  mesh  system) ,  it  is  very  likely  that  the  approxi¬ 
mation  deteriorates  from  third  order  accuracy  in  the  mesh  point 
spacing  to  first  order  accuracy.  In  IV. 18,  the  viscous  stress 
is  expressed  in  terms  of  the  equivalent  rate  of  strain  through 
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equation  IV. 6.  Also  the  expressions  for  the  gradient  of  viscous 
stresses  appear  to  treat  the  dynamic  viscosity  as  independent  of 
x  and  y.  This  is  not  the  case.  Rather  it  can  be  seen  upon 
differentiation  of  the  products  (y(x,y)  for  example,  that 

for  a  power  law  fluid  the  term  (|^)  is  of  higher  order  and, 
therefore,  is  neglected.  Finally  all  terms  in  IV. 19  are  elastic 
xx  stresses.  The  subscripts  and  superscripts  have  been  dropped 
so  as  to  not  severely  encumber  the  equations.  Equation  IV. 18 
is  a  first  order  derivative  counterpart  to  the  steady,  convection- 
dissipation  finite  difference  equation  which  gives  rise  to 
classic  under  and  over-relaxation  methods.  However,  we  do  not 
have  an  equivalent  Courant  number  so  we  merely  employ  Richardson/ 
Jacobi  iteration.  Calling  the  left-hand  side  of  IV. 18  iteration 
k+1  and  the  elastic  stress  terms  on  the  right-hand  side  iteration 
k  (which  is  known)  we  sweep  through  the  entire  solution  domain 
in  the  relaxation  process.  As  in  most  cases,  oe  at  the  first 
iteration  k=l  is  assumed  to  be  zero.  The  issues  which  we  must 
discuss  in  solving  IV. 18  by  this  technique  are  the  selection  of 
mesh  points  i,j,  evaluation  of  the  second  derivative  of  the 
velocity,  convergence  of  the  iteration,  and  treatment  of  boundary 
elements  where  boundary  conditions  must  be  invoked.  We  will  take 
these  in  the  listed  order. 


Since  all  the  terms  involving  the  field  variable  u  in 
equation  IV. 18  are  routinely  calculated,  in  the  evaluation  of 
the  integrals  of  the  finite  element  equations,  at  the  Gauss 
points  in  the  Gauss  quadrature  it  is  natural  to  select  these 
points  as  the  mesh  for  the  elastic  stresses.  Then  for  the 
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differences  required  in  evaluating  the  elastic  stress  gradients, 
the  elastic  stress  at  Gauss  points  of  adjacent  elements  can  be 
used.  This  procedure  is  shown  in  Figure  4  for  one  of  the  Gauss 
points.  Of  course,  two  concerns  arise.  Procedurally,  most 
finite  element  routines  calculate  element  quantities  such  as 
velocity  gradients  in  subroutines  which  dump  the  values  upon 
exiting  the  subroutine,  returning  only  values  of  global  tangent 
stiffness  components.  Therefore,  special  schemes  must  be  devised 
to  identify,  maintain,  and  pass  current  values  of  elastic  stresses, 
external  to  the  subject  element,  to  the  element  for  an  update  of 
the  elastic  stress  at  its  Gauss  points.  Second,  the  discontinuity 
of  stresses  between  elements  which  gives  rise  to  the  practice  of 
"smoothing"  must  be  recognized.  At  the  early  stages  of  itera¬ 
tion,  this  might  aggravate  the  numerical  stability.  For  this 
study,  the  first  issue  was  resolved  by  programming  techniques 
(principly  by  creating  arrays  which  were  stored  in  common 
memories  between  subroutines) .  The  second  issue  was  not  ad¬ 
dressed. 

For  the  problem  of  the  evaluation  of  the  second  derivative 
(recall  from  Section  II  we  are  using  a  "weak"  form  of  the  equa¬ 
tions  so  that  only  CQ  continuity  is  required  of  u) ,  we  now 
require  continuity  of  the  trial  functions  and  explicitly 
evaluate  the  term  just  as  is  done  for  the  first  derivative. 

To  do  this,  a  subroutine  was  written  (ESHAP)  which  returns 
2  ?  2 

the  values  — ?i,  £-Si,  |~ri,  at  the  Gauss  points  of  an  element. 

3x2  3y2  3x37 


The  values  — etc.  are  then  calculated  in  the  exact  same 
3x^ 

manner  as  done  for  the  first  derivatives.  For  this  subroutine 
of  course,  it  was  also  necessary  to  calculate  the  determinant 
of  the  Jacobian  of  the  second  derivatives.  The  mathematics 
involved  in  subroutine  ESHAP  are  given  in  Appendix  2. 

Considering  for  the  time  being  only  convergence  of  the 
Richardson/Jacobi  iteration  scheme,  (Newton-Raphson  and  Picard 
iteration  are  briefly  treated  later) .  We  can  apply  the  Lax/ 
Richtmyer  amplification  matrix  error  method  [31]  discussed  in 
[2].  Briefly,  we  write  equation  IV. 13  in  terms  of  the  final 
value  and  errors  at  each  iteration  or 


(oe  +  e)k+1  =  X  £a(oV  +  oe  +  e)-(u-V) (ov  + 
Subtracting  IV. 13  from  IV. 20  we  get 


(IV. 21) 

(IV. 22) 


ek+1  =  X (A- (u • VT ) £k 


or 


.k+1 


=  X (A  -  u* V) <1 


cre  + 


(IV. 20) 


The  test  for  convergence  then  is  for  the  eigen  values  of  the 
matrix  X (&  -  u*7)  to  be  <1.  Note  that  the  dimensions  of  this 
tridiagonal  matrix  are  3np  where  n  is  the  number  of  Gauss  points 
per  element  and  p  is  the  number  of  elements.  The  complete 
matrix  is  formed  by  assembling  the  individual  3x3  matrices 
at  each  Gauss  point.  We  did  not  conduct  any  further  analysis 
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of  convergence,  but  rather  have  established  bounds  emperically. 
Little  emphasis  was  placed  on  this  issue  because  it  was  found 
during  the  course  of  the  study  that  the  outer  iteration  of 
equation  IV. 12  generally  controlled  convergence. 

Finally  at  the  boundary  elements  where  an  adjacent 

element  may  not  exist,  it  is  necessary  to  devise  an  auxiliary 

3a®  3a® 

scheme  for  the  calculation  of  and  at  the  Gauss  points. 

If  ae  is  known  at  the  element  edges  (in  particular  the  node 
points)  the  nodes  can  be  used  as  the  forward  or  backward  mesh 
points  and  the  relaxation  procedure  continued.  However,  there 
are  some  major  drawbacks  to  this.  First  regardless  of  the 
boundary  condition  (velocity  or  traction  specified)  additional 
calculations  for  velocity  gradients  and  viscous  stress  gradients 
at  the  nodes  must  be  accomplished.  Additionally,  the  elastic 
stress  gradient  can  not  employ  central  differences  at  the  node, 
but  must  be  based  on  a  backward  difference.  Third,  the  forma¬ 
tion  of  the  two  independent  equations  to  simultaneously  solve 
3ae  3ae 

and  is  quite  cumbersome.  A  different  technique  was 
therefore  developed. 

A  new  common  array  was  established  (BOSIG)  to  identify 
and  pass  the  elastic  stresses  at  the  four  corner  nodes.  At 
the  first  iteration,  these  stresses  (four  nodes  by  three  stress 
components  by  the  number  of  elements)  are  initialized  at  zero. 
The  velocity  vector  u  is  then  calculated  in  the  Picard 
iteration.  Then  during  the  calculation  of  element  values  at 
the  Gauss  points  (velocity,  velocity  gradient,  stress  gradients, 
etc.),  the  boundary  elastic  stress  at  the  corner  node  which 
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matches  the  Gauss  point  is  calculated  according  to: 


N.P.  G.P.  .  -e 
— e  -  £e  +  I3T 


G.P. 


,  N.P.  G.P.,  ,  3-e 

(x  "  x  }  +  Ty~ 


(yN-P-.  yG.P. 


G.P. 


(IV. 23) 

Where  N.P.  is  the  node  point  and  G.P.  is  the  Gauss  point. 
This  value  of  elastic  stress  is  then  used  in  the  central  differ¬ 
ence  calculation  at  the  Gauss  point  if  the  element  is  on  a 
boundary.  Figure  5  shows  the  details  for  the  calculation  at  the 
Gauss  points  for  both  boundary  and  interior  elements  as  described 
above. 

To  keep  our  thoughts  clear,  it  is  instructive  to  pause 
and  review.  The  creeping  flow  finite  element  equation  to  be 
solved  is: 


(BT  D  B  +  (mT=1)Ta  mTB)dv  j  u  +J  NTLT£edv  =  J  NTtdA 

v  A 

A 

The  coefficients  of  u  are  linear  and  ae  is  solved  by  successive 

substitution  for  each  value  of  u.  Notice  two  things.  First, 

NTLT  =  BT  so  that  we  could  make  this  substitution.  This  study, 

however,  included  the  terms  Voe  in  the  equation  and  so  these 

m 

values  were  used  directly  with  N  in  calculating  the  integral . 
Second,  a  nested  iteration  on  oe  is  really  not  necessary. 

A 

Rather  we  could  calculate  a  new  u  for  each  update  of  a®  and 
combine  the  two  iterations.  Figure  6  shows  the  two  different 
schemes.  While  not  mathematically  demonstrated,  it  was  felt 
that  such  a  scheme  would  further  degrade  convergence  since 
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u  would  undergo  much  larger  variations.  This  issue  should  be 
considered  in  much  more  depth  in  continuing  studies.  This 
section  will  be  concluded  with  a  discussion  of  three  topics, 
two  very  important,  one  included  only  for  completeness.  These 
topics  are:  convergence  of  the  solutions,  simplication  due  to 
ignoring  the  stress  gradient  terms  of  the  constitutive  equation, 
and  equations  used  for  independently  interpolating  the  total 
deviatoric  stresses  in  a  mixed  finite  element  method.  We 
will  discuss  these  topics  in  order. 

Engelman  et.  al.  [32]  consider  the  problem  of  convergence 
of  the  general  Navier-Stokes  equation  noting  that  Picard  itera¬ 
tion  converges  more  slowly  than  Newton-Raphson ,  but  normally 
over  a  larger  radius.  They  then  treat  the  issue  of  accelerating 
convergence  by  employing  guasi-Newton  methods  emphasizing 
Broyden-Fletcher-Goldfarb-Shanno  updating.  Such  acceleration 
methods  would  enlarge  the  number  of  elements  which  cam  be 
economically  treated  in  the  solution  scheme.  Currently,  however, 
this  is  not  the  problem  with  viscoelastic  flow.  As  we  will  dis¬ 
cuss  in  Section  VI,  the  radius  of  convergence  is  the  major  issue, 
not  the  rate  of  convergence*  Our  study  succeeded  in  obtaining 
solutions  for  Ws^O.Ol  which  could  possibly  be  considered  a  trivial 
case.  However,  for  the  general  flow  geometries,  we  treated  (in 
particular  entry  flow) ,  the  studies  cited  in  the  beginning  of  this 
section  failed  to  achieve  solutions  even  at  that  limit.  Con¬ 
vergence  therefore  is  the  critical  barrier  to  obtaining  more 
general  viscoelastic  solutions.  We  did  not  pursue  such  extensions 


in  this  study,  but  it  is  worthwhile  to  suggest  a  possible 
approach.  Chung's  [2]  review  of  standard  solution  techniques 
is  directly  to  this  point.  The  radius  of  convergence  can  be 
widened  by  continuation  methods.  In  particular,  Chung  suggests 
a  multiple  solution  technique  which  combines  incremental  load¬ 
ing  with  Newton- Raph son  corrections.  Future  effort  in  this 
field  should  investigate  such  an  approach.  We  employed  Picard 
iteration  exclusively.  Picard  iteration  should  be  tried  as  the 
top  level,  along  with  continuation  methods.  It  is  noted  that 
both  types  of  solution  are  amenable  to  the  computer  program 
used  in  this  study. 

We  turn  now  to  the  simplications  when  the  stress  gradient 
terms  are  neglected.  The  terms  themselves  arise  in  the  con¬ 
vection  terms  of  the  constitutive  equation,  i.e.,  (u*V)o.  For 
creeping  flow  similar  terms  were  neglected  in  the  Navier-Stokes 
equation  and  we  know  that  for  polymer  melts,  this  is  a  good 
approximation.  It  is  then  obvious  that  we  compare  approximate 
magnitudes  of  Vu  and  Vo.  For  viscoelastic  flows,  we  have  already 

©  V 

established  that  o  is  on  the  order  of  a  and  the  gradients 
might  be  expected  to  be  of  equivalent  nature.  Therefore,  we 
look  at  the  comparison  between  the  first  derivative  of  u  and 
the  second  derivative.  It  is  known  that  even  when  u  is  discon¬ 
tinuous  (as  in  the  case  of  cross-channel  flow  of  a  screw  extruder 
[9],  the  approximation  at  small  distances  from  the  singularities 
of  Vu  are  quite  good.  This  suggests  that  for  creeping  flow,  a 
good  approximation  may  be  achieved  when  (u*V)c  is  neglected. 


Equation  IV. 13  then  becomes: 

£*  -  XA(av  +  ae)  (IV. 24) 

Evaluation  of  7a  is  eliminated  and  the  Picard  iteration  be¬ 
comes  much  more  straightforward.  An  optional  approach  is  to 
solve  ae  explicitly  as: 

(I  -  XA)oe  ■  XA  aV  (IV. 25) 

or 

£6  *  U  -  XA)'1  XA  aV  (IV. 26) 

Where  I  is  the  unit  matrix  6 i j  ^j1  (i,j  =  1,2,3) 

Equation  IV. 26  allows  IV. 5  to  be  written  explicitly  in  terms 

a 

of  u  and  the  equation  is  a  simple  non-linear  equation  which 
can  be  solved  with  the  numerical  techniques  discussed  throughout 
the  report.  It  is  noted  that  although  the  explicit  form  makes 
the  equations  more  straightforward,  it  is  not  expected  that 
the  radius  of  convergence  (which  is  a  function  of  X)  will  be 
widened  much.  However,  at  the  early  stages  of  research  efforts, 
particularly  in  applying  continuation  methods,  this  equation 
seems  to  offer  promise. 

Finally,  the  mixed  method  of  solution  is  briefly  discussed 
for  sake  of  completeness.  Following  Kawahara's  approach  [28], 
we  set  up  the  simultaneous  equations  for  steady  state  in  indicial 
notation: 
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p  u.  uifj  +  p,i  -  OLj ,  j  = 


(IV. 27a) 


°ij  +  X(uk0ij,k  -  ui,k0kj  '  "j^ik1  -  +  uj#i>  “  ° 

(IV. 27b) 

Both  IV. 27a  and  IV. 27b  are  non-linear;  we  write  the  finite 


element  equations : 


^(NTP(V*(N  u)T)TN  +  (mT  B)TamT  B)dv 


f  *fp  *rp 

/( N  iXV(N  a)N  -  N  D  B  -  N  Q  N)dv 


u+  |yNTB*Tc 

u  +|yN*TN*c 
**v 


(IV. 28a) 

VTav]  ;  -  t 


dv  a  =  o 


(IV. 28b) 


Where  the  asterisk  indicates  the  trial  function  for  the 
stress  interpolation. 

The  solution  to  IV. 28  can  be  seen  clearly  if  we  form  a  typical 
equation  in  matrix  form: 

NTp  (V*  (N  u)  )  N.  i  UA  F 

=  —  “3  i  i 

i 

T  T  T  I  T  *T  *1  —1 

+  <E  Si>  »E  fij  !  Ni  Bj  v  F2  (IV.  29) 


*>p  */v 

XV  (N  o)N- 

*T 

-Ni  2  lj 

-n*T£  n. 

i  *  =*i 


(IV. 29) 


*»r  * 
Ni  Nj 


(In  equations  IV. 29,  the  integrals  are  implied.) 
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In  IV. 29,  the  superscript  in  the  column  vectors  indicate  the 
node  number  so  that  this  relation  is  repeated  for  each  of 
the  nine  nodes,  i  and  j  indicate  the  row  and  column  in  the 
assembled  array  (for  IV. 29  i=j=l) .  The  array  is  partitioned 
accordingly  so  that  the  upper  left  corner  is  2  x  2,  upper  right 
corner  is  2  x  3,  lower  left  corner  is  3x2,  lower  right  corner 
is  3  x  3.  All  matrices  in  IV. 29  have  been  previously  defined 
with  the  exception  of  Q  which  is: 

2(s  “xx>55  +  2(S  «xy>!y 

2=>  o  2<S*Vfe+2 

a  */v  a  a  * 

NO  £-  +  N  a  4-  No  4—  +  N 

-  yy  3y  -  xy  9x  -  xx  8x  =- 

These  equations  when  fully  assembled  yield  a  set  of  linear 
equations  of  order  5p,  where  p  is  the  number  of  nodes.  For  a 

nine  node  element  then  the  order  of  equations  is  45.  The  number 
of  variables  for  the  whole  domain  then  would  be  45n-m  with  n 
being  the  number  of  elements  and  m  the  number  of  shared  nodes. 

It  can  be  seen  that  it  does  not  take  many  elements  to  generate 
a  very  large  computer  region  to  solve  the  equations.  While  the 
above  analysis  was  conducted  and  subroutine  EIMT06  written  for 
the  problem  solution,  no  flow  cases  were  run  in  this  study. 
Future  work  may  implement  subroutine  ELMTJJ6. 


V.  COMPUTER  IMPLEMENTATION 


In  this  section  we  will  discuss  the  major  aspects  of  the 
finite  element  program,  the  calculational  procedures,  and  the 
input/output. 

The  source  program  was  a  modified  version  of  the  Finite 
Element  Analysis  Program  (FEAP)  written  by  Prof.  R.  L.  Taylor 
at  the  University  of  California,  Berkeley,  and  published  in 
Chapter  24  of  [l].  The  modifications  have  been  made  by 
Prof.  David  Roylance  of  the  Massachusetts  Institute  of  Tech¬ 
nology  to  accommodate  polymer  melt  flow  [9].  These  modifica¬ 
tions  are  largely:  (i)  addition  of  a  power  law  flow  rule, 

(ii)  addition  of  a  temperature  dependent  viscosity,  (iii) 
alteration  of  matrix  algebra  operations,  and  (iv)  addition 
of  an  axisymmetric  capability.  The  rationale  for  using  this 
model  is  given  in  [9].  The  current  effort  included  reviewing 
the  source  program  to  insure  correctness,  and  modifying  it 
to  include  a  viscoelastic  flow  option.  Currently  the  program 
is  two-dimensional  (rectilinear  or  axisymmetric)  and  steady 
state . 

The  program  establishes  a  dynamic  storage  vector  at  the 
outset  which  is  partitioned  to  store  all  input  data  (node  co¬ 
ordinates,  element  node  numbers,  etc.),  global  data  (stiffnesses, 
loads,  etc.)  and  output  data  (velocities).  Other  features  are 
a  linear  interpolation  mesh  generation  scheme,  an  active 
equation  solver  and  a  macro  command  language  which  controls 
the  solution  execution.  The  macro  commands  and  their  meaning 


are  listed  in  Table  24.12  of  [l]. 

Upon  construction  of  the  architecture  of  the  problem, 
calculations  required  for  a  specific  command  (such  as  forming 
the  tangent  stiffness  matrix)  are  made  in  a  library  of  element 
subroutines.  Subroutine  PFORM  steps  through  the  n  elements  by 
forming  element  arrays  from  global  data  and  passing  the  arrays 
to  the  element  routine.  Subroutine  ELMT05  is  a  general  2D 
penalty  method  solution  of  the  Navier-Stokes  equation  written 
by  Frecaut  [16].  This  is  the  element  subroutine  modified  for 
the  viscoelastic  flow. 

The  basic  source  program  flow  chart  is  given  in  Figure  7 . 

To  modify  this  program  for  viscoelastic  flow,  three  basic 
changes  were  made.  First  was  to  flag  the  problem  as  visco¬ 
elastic  and  read  material  data.  This  was  done  in  subroutine 
DFMTRX.  The  card  reading  format  after  input  macro  command 
MATE  was  changed  to  the  following: 

CARD  1  Format  (15,  4X,  II,  17A4) 

CARD  2  Format  (415,  F10.0) 

CARD  3  Format  (15,  7D1JJ.4) 

Card  one  reads  the  material  set  number  in  columns  1-5  (in  all 
cases  only  one  material  set  is  used  and  therefore  this  is  1) , 
the  element  type  in  column  10  (5  for  ELMT05)  and  the  problem 
description  in  the  remaining  coiumns.  Card  two  reads  the  flow 
type  in  columns  1-5  (1  =  plane  flow,  3  =  axi«?ymmetric  flow)  , 
a  flag  (Nl)  for  thermal  coupling  in  columns  6-10  ( 0  =  isothermal, 
1  =  thermally  coupled) ,  a  flag  (K2)  for  viscoelasticity  in 
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columns  11-15  (1  =  simple  viscous,  2  =  power  law  viscous, 

3  =  White-Metzner  Viscoelastic,  4  =  DeWitt  Viscoelastic, 

5  =  Rivlin  Ericksen  Viscoelastic) ,  a  flag  (N3)  for  the  time 
domain  in  columns  16-20  (1  =  steady  state,  2  =  unsteady), 
and  the  power  law  coefficient  (P4)  in  columns  21-30.  P4 
must  be  included  and  for  simple  viscous  material  P4  =  1.0 
(which  was  the  case  treated  exclusively  in  this  study) . 

Card  three  reads  the  Gauss  integration  order  (L)  in 
columns  1-5  (2  =  2x2) ,  the  penalty  coefficient  (XLAM)  in 
columns  6-15,  the  viscosity  coefficient  (XMU)  in  columns 
16-25,  the  density  (RHO)  in  columns  26-35,  the  viscoelastic 
shear  modulus  (G)  in  columns  36-45,  the  thermal  conductivity 
(XK)  in  columns  46-55,  the  specific  heat  (C)  in  columns 
56-65,  and  the  work-to-heat  conversion  factor  (HEAT)  in 
columns  66-75.  The  program  is  written  so  that  when  data  is 
not  required  for  the  specific  problem  (e.g.  linear,  steady, 
isothermal,  inelastic  flow)  those  columns  may  be  left  blank. 

In  card  three  then  only  columns  1-25  need  be  included. 

The  second  change  was  to  add  algorithms  in  ELMT05  for 
the  calculation  of  the  elastic  stresses  according  to  equations 
IV. 13.  The  last  change  presented  the  major  difficulty:  the 
calculation  of  the  elastic  stress  gradients  according  to 
equations  IV. 19.  As  noted  in  the  previous  section,  no  scheme 
existed  for  making  calculations  with  variables  from  different 
elements.  In  order  to  solve  IV. 19,  however,  this  was  necessary. 

The  approach  taken  was  to  define  common  arrays  YY(I,J,N),  ESIG1(I,J,N) 
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ES1G2 (I , J, N) , ESIG3 (I, J,N) , ELAS1 (I , J,N) ,ELAS2 (I , J,N) , ELAS3 ( I , J,N) , 

and  BOSIG(I,J,N) .  YY  is  the  global  coordinate  (J=l,2)  of  the 

Gauss  points  (1=1,4).  ESIG1,  ESIG2 ,  and  ESIG3  are  a  ,  a 

xx  yy 

and  a  respectively  at  the  Gauss  points  (1=1,4)  at  iteration 

Ay 

J=K,  K+l.  ELAS1,  ELAS2 ,  and  ELAS3  are  the  gradients  (J=l,2) 
of  oe  at  the  Gauss  points  (1=1,4) .  BOSIG  is  the  elastic  stress 
(J=l,3)  at  the  boundary,  at  the  corner  node  (1=1,4).  In  all 
the  arrays  N  is  the  element  number.  In  PFORM,  N  is  passed  as 
common  through  ELMLIB  and  ELMTJJ5  and  it  is  therefore  possible 
to  conduct  the  calculations  between  the  two  subroutines  PFORM 
and  ELMT05.  The  gradients  of  the  three  stress  components  at 
the  Gauss  points  are  first  solved  for  all  the  elements  assuming 
they  are  a  boundary  element  on  all  sides.  A  searching  scheme 
is  then  affected  which  compares  the  nodes  of  all  the  other 
elements.  When  two  elements  are  found  in  the  correct  location, 
the  elastic  stress  gradients  are  replaced  at  that  Gauss  point. 

If  adjacent  elements  are  not  found,  the  element  is  on  a  boun¬ 
dary  and  that  Gauss  point  is  left  unchanged.  During  the 
Richardson/ Jacobi  iteration,  the  elastic  stress  gradients  then 
are  calculated  in  PFORM  and  these  values  used  in  ELMT05  to 
calculate  the  updated  values  of  the  elastic  stress  at  the 
K+l  iteration.  This  iteration  is  conducted  20  times  unless 
convergence  is  achieved  beforehand.  The  program  then  continues 
in  a  normal  manner. 

The  listings  of  the  major  subroutines  written  to  accomplish 
the  modifications  are  included  in  Appendix  3.  The  subroutines 
are  in  order  listed  ELMT05,  ELMT06,  ESHAP,  PFORM,  CMATRX,  and 
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FPSIG.  ELMT06  is  the  subroutine  written  for  interpolate  total 
deviatoric  stresses  in  a  mixed  method.  ESHAP  is  the  calcula¬ 
tion  of  the  second  derivatives  and  FPSIG  is  a  new  routine 
written  to  print  viscous  and  elastic  stresses  at  the  Gauss 
points.  CMATRX  is  the  subroutine  which  forms  the  Q  matrix 
in  ELMT06. 
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VI .  CALCULATION  RESULTS 


Four  flow  geometries  were  treated  as  shown  in  Figure  8 
(along  with  the  boundary  conditions) :  Cross  Channel  Flow, 

Plane  Couette  Flow,  Entry  Flow,  and  Step  Flow.  Table  1  shows 
the  computer  run  matrix.  The  input  data  sets  for  runs  1,  3,  4, 
6,  13,  and  20  are  included  as  Appendix  4.  Results  are  dis¬ 
cussed  below  for  each  of  the  four  problems  treated.  For  all 
cases,  the  viscosity  coefficient  was  taken  to  be  790  poise. 

This  was  the  value  selected  by  Roylance  [9]  in  previous  studies. 
His  reasons  were  unrelated  to  the  work  in  this  study,  but  we 
chose  to  use  the  same  value  for  comparison  purposes.  With 

4 

more  reasonable  values  (10  ) ,  we  would  only  expect  to  see 
higher  stresses,  but  no  change  in  the  velocity  fields. 

CROSS  CHANNEL  FLOW 

The  solution  of  creeping  flow,  circulating  in  the  trans¬ 
verse  plane  of  channel,  for  a  viscous  fluid  is  well  known  (e.g. 
[9]).  At  steady  state,  the  circulation  is  uniform  with  a  vortex 
center  at  mid-height,  towards  the  vertical  boundary  on  the  right 
in  Figure  8a.  This  study  looked  at  the  consistency  of  repro¬ 
ducing  this  flow  with  9  node  and  8  node  elements  and  the  effects 
of  a  finite  fluid  elasticity.  Secondary  eddies  and  screw  power 
requirement  changes  were  considered  to  be  demonstrable  effects 
of  elasticity. 
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Figure  9  shows  the  velocity  vector  flow  field  for  run  1 
(linear  case).  Results  are  identical  to  [7],  different  from  [9]. 
This  is  due  exclusively  to  the  specified  boundary  condition  at 
the  upper  corners  of  the  channel.  For  our  boundary  conditions, 
the  vortex  center  is  at  the  mid-width  of  the  channel  near  the 
2/3  height  section. 

The  velocities  calculated  for  the  nodes  of  elements  7,  9, 
and  15  by  the  18  element  9  node  and  18  element  8  node  case 
are  compared  in  Figure  10.  Note  that  a  significant  difference 
occurs  in  the  direction  of  the  resultant  velocities  in  element  7 
and  the  magnitude  in  element  15  (a  20%  lower  horizontal  velocity 
is  predicted  in  the  middle  nodes  of  element  15  by  the  8  node 
model) .  When  the  results  of  the  72  element,  8  node  case  are 
examined  (run  3)  the  9  node  model  is  found  to  be  uniformly 
closer.  The  velocity  field  is,  therefore,  predicted  much  better 
by  the  9  node  elements  for  the  same  number  of  elements. 

Let  us  now  make  a  practical  application.  The  power  per 

unit  area  required  of  a  single  flight  screw  extruder  to  create 

this  circulation  is  the  shear  stress  in  the  fluid  times  Ufi 

(the  relative  barral  velocity) .  If  we  approximate  this  as  the 

average  element  shear  stress  a  times  the  average  velocity 

xy 

in  the  element,  we  have  the  following  for  element  15: 


9  NODE  18  ELEM  8  NODE  18  ELEM 


O  (dynes/cm2)  0.22  x  106  0.2  x  106 

xy 

u  (cm/sec)  -50  -52.5 

w  (dyne-cm/cm2-sec)  1.08  x  10^  1.05  x  10^ 

We  can  conclude  that  the  9  node  elements  yield  more  accurate 
node  velocities,  but  when  average  properties  are  sought,  such 
as  the  power  or  torque  required  for  the  screw  design,  both 
models  give  approximately  the  same  results  for  equivalent 
meshes.  This,  of  course,  is  expected  since  the  finite  element 
equations  satisfy  equilibrium  over  the  entire  region.  However, 
on  a  local  scale  (which  we  are  also  interested  in)  the  above 
justifies  our  earlier  preference  for  the  9  node  elements. 

From  Hughes  data  [7],  the  effects  of  increasing  the 
Reynold’s  number  (Re)  is  to  shift  the  vortex  center  toward  the 
right-hand  boundary.  This  was  investigated  for  one  case  by 
choosing  the  density  of  polyphenelenesulfide  (1.6  gm/cm3) . 
Combining  this  with  the  other  characteristic  numbers  of  the 
cross-channel  flow  problem,  we  obtain  Re  =  =  0.41. 

Including  the  convection  non-linearity  for  this  Re  we  found  no 
discernible  perturbation  to  the  velocities  or  stresses,  thus 
confirming  the  validity  of  the  "creeping  flow"  analysis. 

For  the  single  viscoelastic  case  for  which  the  solution 
converged  (Ws  =  0.02)  the  velocity  field  again  did  not  vary 
appreciably.  Figure  9  can,  therefore,  be  considered  correct 
for  this  level  of  elasticity.  To  look  at  the  stress  effects. 


we  make  the  same  calculation  for  the  specific  power  as  above 
yielding: 

NEWTONIAN  VISCOELASTIC  (W s-0.02) 

o  (dynes/cm2)  0.22  x  106  0.22  x  106 

xy 

u  (cm/sec)  -50  -50 

*  ?  7  7 

w  (dyne-cm/cm  -sec)  1.08  x  10  1.08  x  10 

Within  roundoff  error,  the  two  flows  are  identical  (maximum 

a  deviation  was  1%) .  A  second  comparison  is  available  in 
xy 

Figure  11  where  the  pressure  is  plotted  at  the  mid-height  as 
a  function  of  the  cross-channel  (transverse)  station.  Again 
the  viscoelastic  flow  is  coincident  with  the  Newtonian  case. 
Within  the  range  of  calculations  achieved  in  this  study 
therefore  (Ws<0.02)  ,  there  are  no  effects  of  viscoelasticity 
manifested.  We  do  observe,  however,  that  the  stresses  calcu¬ 
lated  (-1%  variation)  are  consistent  with  the  Ws  suggesting 
accuracy  of  the  computer  model  when  convergence  is  achieved. 

PLANE  COUETTE  FLOW 

Plane  Couette  flow  was  selected  for  the  fundamental 
evaluation  of  the  computer  model.  This  is  through  the  relation 
presented  by  Middleman  [26]: 

SR  ' 

where  SR  is  the  recoverable  (elastic)  shear  stress: 


(VI. 1) 


X  is  the  relaxation  modulus  and  y  is  the  steady,  simple 
shear  flow  strain  rate.  The  flow  is  enforced  by  specifying  a 
linear  variation  of  the  horizontal  velocity  between  two  plates, 
one  stationary,  the  other  moving  at  a  constant  velocity  as  shown 
in  Figure  7b. 

Run  7  was  the  Newtonian  case  to  validate  the  problem.  In 

this  case,  cr  „  and  o„„  should  be  identically  zero  and 
xx  yy 

constant  throughout  the  field  domain.  This  was  the  result  of 
the  calculation. 


For  the  viscoelastic  case  (Run  10) ,  all  the  normal  stresses 

£ 

are  elastic  while  from  equations  IV. 13,  with  v  =  ^  *  °  only 
is  finite.  Therefore,  we  should  observe  the  following: 


c  _  xx 


-  Xy  =  constant 


(VI. 3) 


For  a  unit  height  between  sliding  plates  we  have  y  =  UB  so 
that: 

a  ®  =  2XUn  o  V  (VI. 4) 

xx  B  xy 

The  computer  results  are  for  X  =  0.0002,  UB  =  100  cm/sec 
(Ws  =  0.02)  : 


o^®  =  3.16  KPa,  2XUBaxyv  =3.16  KPa. 

The  equation  is  identically  satisfied.  This,  of  course,  is 
encouraging  for  future  work  to  increase  the  radius  of  convergence 
for  higher  Ws  numbers . 
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ENTRY  FLOW 


The  entry  flow  problem  for  viscoelastic  fluids  has  not 
been  successfully  calculated  by  finite  element  methods  in 
the  past,  due  to  severe  numerical  convergence  problems.  As  a 
first  step.  Run  11  was  accomplished  for  linear  flow  according 
to  the  boundary  conditions  specified  in  Figure  7c.  A  discus¬ 
sion  of  these  boundary  conditions  is  in  order. 

Rather  than  a  constant  horizontal  velocity  at  the  inlet 
to  the  reservoir  (upstream  channel) ,  a  more  accurate  analysis 
would  specify  fully  developed  flow.  Middleman  [26]  presents 
this  for  flow  between  parallel  plates  (for  a  Newtonian  fluid 
as)  : 

3P 

“  -  frr  [*  -  (s')  ]  (VI-5’ 

where  B  is  the  channel  height 
and  L  is  the  channel  length 

(all  other  variables  retain  their  earlier  definition) . 

For  a  White-Metzner  fluid,  the  plane-Poiseville  flow 

would  be  solved  by  adding  the  elastic  stresses  to  the  momentum 

equations.  Perera  [27]  did  this  for  a  4  constant  Oldroyd  fluid 

and  solved  the  resulting  second-order  differential  equation  for 

u(y)  by  Newton-Cotes  integration.  With  equations  of  the  type 

dP 

specified  in  VI. 5,  we  can  solve  the  pressure  loss  ^  due  to 
inlet  and  outlet.  In  addition  White  [33]  cites  the  additional 
pressure  losses  due  to  entrance  and  exit  of  the  dies.  It  is 
these  boundary  conditions  that  would  be  more  realistic  in 
treating  the  entry  flow  problem  (velocity  according  to  VI. 5 


at  one  end,  at  the  other) .  With  the  formulation  specified 
in  this  work,  it  was  expected  that  the  flow  field  would  behave 
quite  differently  from  the  classical  converging  type.  Since 
we  did  not  have  data  on  die  pressure  losses,  however,  the 
initial  calculations  were  made  on  the  basis  of  the  boundary 
conditions  given. 

When  fully  developed  conditions  are  specified,  both  upstream 
and  downstream  of  the  entrance  region  the  flow  is  known  to  be 
stable  up  to  relatively  high  Ws  numbers .  At  Ws  around  one 
secondary  vortex  patterns  arise  which  are  generally  ascribed 
to  increasing  elastic  stresses  generated  in  the  shearing/elonga- 
tional  flow  (White  [33]  implies  that  elongational  flow  is  im¬ 
portant  and  we,  therefore,  conclude  that  the  Rivlin-Ericksen 
fluid  simplified  for  viscometric  flow  is  a  questionable  model) . 
This  flow  behavior  is  documented  in  Figure  12  which  shows 
experimental  behavior  noted  by  White  [33]  as  a  function  of  Ws 
and  calculations  of  Perera  [27]  for  Ws  =  0.6. 

The  calculated  velocity  field  for  the  boundary  condition 
specified  in  Figure  8c  is  shown  in  Figure  13.  Although  the 
mesh  is  very  coarse,  it  appears  that  the  flow  is  unstable  for 
these  conditions.  The  viscoelastic  calculation  (Run  12,  Ws  = 
0.01)  exhibited  identical  behavior.  Because  of  this  poorly 
behaved  flow  field,  the  calculation  was  repeated  using  the 
fully  developed  flow  boundary  conditions.  The  results  are 
shown  in  Figure  14.  The  specific  boundary  conditions  were 
established  in  the  following  manner.  The  excess  pressure 
losses  described  by  White  [33]  were  ignored  (this  will  affect 
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the  calculation  however) .  At  y  *  0  (y  measured  from  the 
mid-height  of  the  channel)  equation  IV. 5  is 


u  = 


b2ap 

8yL 


(IV. 6) 


For  the  two  channels,  there  would  be  a  total  pressure  loss  of 
APt  *  APj  +  APq  if  the  flow  was  fully  developed.  Therefore 


L  u  Lu 

AP_  =  8y 

BI  Bo 


(IV. 7) 


For  our  geometry  Lj  =  L  =  B^.  =  1,  Bq=  Jj,  and  p  =  790  poise. 
There  are  three  unknowns  in  equation  IV. 7.  However,  rather 
than  specifying  two  of  the  three,  we  merely  let  APj  =  APq 
(which  has  the  same  effect)  and  specified  uq.  This  permits 
the  calculation  of  u^.  and  thereby  calculation  of  APT.  This 
APt  was  established  as  the  inlet  traction  and  the  outlet 
was  atmospheric  PQ  ■  O.  This  yields  the  value  of  Pj  *  6.4p. 

The  pressure  is  converted  to  the  virtual  "work"  equivalent 
node  forces  by  the  relationship 

FC  =  i  H  PT 
x  3  I 

pM  -  i  b  p 

x  31  (IV. 8) 

where  the  superscript  denotes  the  element  node  (c  -  corner 
node,  M  =  mid- side  node)  and  H  is  half  the  element  height  at 
the  nodes.  (See  Frecault  [16]  for  the  details  of  virtual  "work' 
equivalence  calculations;  equation  IV. 8  are  valid  for  8  node 
and  9  node  plane  quadrilateral  elements) .  (In  the  actual 
boundary  node  forces,  the  corner  nodes  are  loaded  with 


c  2 

Fx  "  3  H  pi  for  uniform  meshes  since  the  node  is  shared  by 
adjacent  elements.  Only  the  vertical  velocities  at  the 
boundaries  (v=o)  now  must  be  specified.  The  mid-height  hori¬ 
zontal  velocity  will  not  be  the  value  used  in  the  calculation, 
but  the  flow  will  be  fully  developed. 

Comparing  Figures  13  and  14,  we  see  that  although  the 
behavior  is  somewhat  improved  by  the  fully  developed  flow  case, 
there  is  still  major  error  in  the  flow  field  and  even  flow 
reversal.  This  is  felt  to  be  attributable  entirely  to  the 
coarseness  of  the  mesh,  particularly  near  the  entry  corner. 

A  finer  mesh  case  was  not  constructed  to  test  this  hypothesis. 
It  is  recommended  that  future  work  include  this  refinement. 

Notice  that  symmetry  was  not  employed  to  reduce  the 
number  of  elements.  This  was  due  to  the  difficulty  of  speci¬ 
fying  boundary  conditions  on  the  plane  of  symmetry.  The  first 
conditon  is  v=o,  but  the  other  boundary  condition  is  not  so 

straightforward.  We  know  that  and  are  zero  at  the  mid- 

3  v 

height,  but  in  general  is  not  zero  within  the  reduction 

ay 

region.  This,  of  course,  is  a  statement  that  the  one  dimen¬ 
sional  lubrication  theory  is  not  valid.  Since  the  pressure 
now  changes  across  the  channel  height  the  pressure  in  the  x 
direction  can  no  longer  be  specified  as  a  linear  function  of  x. 
Therefore,  the  nodal  loading  in  the  x  direction  is  unprescribed 
as  well  as  the  velocity  u.  Of  course,  this  could  be  resolved 
by  adding  the  condition  of  no  mass  flow  across  the  plane  of 
symmetry.  We  chose  not  to  accomplish  this  at  the  penalty  of 
doubling  the  number  of  elements. 
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If  the  flow  were  one  dimensional,  the  pressure  would  vary 
linearly  with  the  length  and  the  velocity  would  be  constant 
in  each  of  the  two  sections  (plane  Poiseulle  flow) .  Figure  15 
shows  the  deviation  from  this  case. 

It  was  noted  in  examining  the  stresses  in  Runs  11  and  12 
that  the  difference  was  much  larger  than  expected  for  the  low 
Ws  •  However,  a  thorough  evaluation  was  not  conducted  because 
it  was  felt  that  the  differences  were  an  artifact  of  the 
calculations  due  to  the  following:  (i)  the  velocity  fields 
were  erratic  as  previously  mentioned  due  to  the  coarse  grid, 
(ii)  the  boundary  conditions  of  constant  inlet  velocity  gave 
rise  to  poorly  behaved  pressure  variations  even  for  the 
Newtonian  case,  and  (iii)  the  solution  convergence  for  the 
non-linear  problem  was  still  poor  at  30  iterations.  It  is 
noted  in  passing,  however,  that  as  the  solution  procedure  is 
improved,  it  is  exactly  these  types  of  variations  which  are 
being  sought. 

STEP  FLOW 

This  geometry  was  selected  as  the  beginning  step  toward 
an  analysis  of  flow  past  an  obstruction  such  as  would  be  the 
case  if  pins  were  added  to  the  cavity  to  form  holes  in  the 
molded  part.  With  the  boundary  conditions  specified  in 
Figure  8d,  the  results  were  very  similar  to  those  discussed 
for  entry  flow.  A  discussion  of  the  computer  calculations 
will  therefore  not  be  included  in  this  report.  It  is  noted, 
however,  that  there  is  still  negligible  differences  between 
Runs  15  (linear  Newtonian)  and  16  (convection  Newtonian) . 
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This  begins  to  address  the  issue  of  "Stokes  paradox"  and  the 
necessity  of  including  convection,  even  for  low  Re,  for 
obstructed  flow.  The  paradox  is  that  in  two-dimensional 
flow  no  analytic  solution  exists  for  the  linear  equation 
which  matches  the  boundary  conditions  at  the  surface  of  the 
obstruction  and  at  large  distances  away  from  it  [34], 
Batchelor  [35]  shows  that  when  the  distance  from  the  obstruc¬ 
tion  (or  a  boundary)  is  on  the  order  of  i/Re  (where  t  is  a 
characteristic  dimension  of  the  obstruction)  the  convection 
stresses  (inertia)  may  become  of  equal  importance  to  the 
viscous  stresses.  Analytically  this  correction  is  known  as 
O' seen' s  improvement.  Again  as  the  model  described  in  this 
report  is  refined,  the  adequacy  of  the  "creeping"  flow 
analysis  must  be  examined  in  light  of  this  issue. 
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VII.  MODEL  EVALUATION 

It  is  worthwhile  to  complete  a  qualitative  evaluation 
of  the  computer  model  before  this  report  is  concluded. 

Figure  16  presents  a  diagram  of  a  complete  model  for  a  real 
injection  molding  process.  The  Figure  emphasizes  those 
elements  included  in  this  study.  Since  we  achieved  numerical 
convergence  for  Ws  £  0.01  it  must  be  concluded  that  a  non- 
Newtonian  power  law  fluid  analysis  would  be  as  good  an 
approximation  as  the  viscoelastic  model  used  herein.  If 
future  work  does  not  improve  this  convergence  region  (at 
least  to  Ws  0.5)  the  numerical  analysis  would  seem  to  be 
as  good  without  including  elastic  effects.  Also  finite 
difference  methods  have  succeeded  in  obtaining  solutions  up 
to  Ws  -  0.6  [27]  and  it  may  therefore  be  advisable  to  develop 
these  techniques  for  application  to  the  gyroscope  manufacturing. 

The  model  is  steady  state  and  includes  no  free  surfaces 
such  as  would  occur  during  the  mold  filling  period.  Therefore, 
it  can  only  be  used  in  regions  such  as  the  extruder,  nozzle, 
sprue,  runner  and  gate.  Unless  unsteady,  free  surface  terms 
are  included,  this  model  is  not  applicable  to  the  mold  filling 
itself.  But  the  power  required  to  supply  a  nozzle  with  a  given 
rate  of  flow  is  certainly  within  the  capability  of  the  model. 
Also  the  state  of  the  bulk  material  as  it  passes  through  the 
gate  can  be  determined  by  use  of  this  model.  Any  damage  due 
to  high  stresses  or  thermal  degradation  in  these  regions  can 
be  analyzed  with  the  model.  It  is  noted  that  although  there 
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will  be  a  finite  elastic  stress  which  the  polymer  can  sustain 
before  flowing  completely  plastic  (viscous  plus  the  elastic 
limit  stress),  there  is  no  yield  stress  built  into  this  model. 
Therefore,  while  the  model  will  predict  continually  increasing 
stresses,  judgement  must  be  exercised  as  to  the  real  elastic 
capacity  of  the  fluid. 

The  current  status  of  the  coupled  heat  transfer  capability 
of  the  model  is  the  adiabatic  model  developed  by  Roy lance  [9]. 
Extension  to  a  complete  non-isothermal  boundary  analysis  can 
be  implemented  without  too  much  difficulty. 

We  have  noted  that  major  modifications  are  necessary  to 
evaluate  the  mold  filling  itself  (only  pressure  and  filling 
rate  can  currently  be  analyzed) .  Also  within  the  mold,  the 
cooling  stage  of  the  molding  process  can  not  be  analyzed 
because  of  the  absence  of  a  solid  thermomechanical  viscoelastic 
model. 

However,  if  an  initial  state  can  be  established  for  the 
cooling  process  such  a  model  could  be  developed. 

The  mold  filling  process  itself  can  take  the  approach  of 
a  constant  flow  rate  at  the  gate  once  free  surface  effects  are 
added  to  the  model.  This  is  the  approach  used  in  [11].  The 
free  surface  analysis  is  most  clearly  discussed  in  [4]  where 
the  front  displacement  is  calculated  over  some  interval  of 
time  assuming  a  constant  velocity  of  the  boundary  elements 
node  points.  The  surface  traction  on  the  flow  front  is  zero 
normal  to  the  surface  and  the  material  surface  tension  tangen¬ 
tial  to  the  surface. 
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From  Section  V,  we  can  discuss  the  approach  to  improv¬ 
ing  the  viscoelastic  case.  To  assess  the  maximum  radius  of 
convergence  of  the  momentum  equation,  it  is  adequate  to 
neglect  u*Vc  and  use  equation  IV. 26.  Since  Newton-Raphson 
iteration  generally  converges  for  the  Navier-Stokes  equations 
well  above  Re  =  25,  it  should  be  verified  that  convergence 
is  achieved  with  the  current  numerical  approach  for  Ws  =  25. 
With  this  step  accomplished,  u*Vo  can  be  added  and  the  con¬ 
tinuation  method  used.  The  effective  technique  should  employ 
incremental  loading  with  Newton-Raphson  corrections.  Let 
us  discuss  this  a  little  further.  Since  we  are  using  direct 
(Picard)  iteration  on  the  elastic  stress  terms,  let  us  rewrite 
equation  IV. 11  as: 

K  u  =  f  -  K0  (VII. 1) 

Since  Picard  iteration  is  a  single  point  scheme,  (i.e., 

A 

the  initial  value  of  K  u  +  if-  f  is  always  used  rather  than 
updating  in  the  Newton-Raphson  scheme,  see  Figure  17)  we  can 
attempt  to  increment  this  point.  Therefore,  instead  of  solving 
VII. 1  directly,  we  solve: 

K  u  =  9(f  -  if  )  (VII. 2) 

where  0<Q<1.  With  the  solution  to  VII. 2  converging  for 
sufficiently  small  numbers  of  0  we  can  update  the  initial  se- 

A 

lection  of  u  by  incrementing  0.  For  example,  let  0=0.1  then 

/v 

in  the  first  increment  the  first  value  of  u  is  : 
u°  =  K-1  0f 
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We  then  iterate  with  K  us+1  =  e[F  -  K® (us) ] . 

When  convergence  is  achieved  then  we  increment  8  to  2©  =  0.2. 

Then  ~  ,  _  „ , . 

u  =  K_1[2©F  -  K  (u+1) ] 

Therefore,  the  initial  guess  is  improved  by  the  correction 
0  s+1 

K  (u  ) .  It  is  noted  that  this  technique  is  different  from 
the  normal  continuation  methods  where  the  non-linear  equation 
is  always  of  the  form:  K(u)u=f.  While  no  mathematical  analy¬ 
sis  has  been  conducted  on  this  proposed  technique,  it  appears 
to  offer  promise. 

This  deviation  in  the  classical  incremental  load  method 
is  only  necessary  when  the  stress  gradient  terms  are  included 
in  the  viscoelastic  constitutive  model.  Therefore,  when  the 
model  undergoes  its  first  revision  with  u*Va  neglected  we 
write  oe  explicitly  and  if  convergence  fails  the  classical 
incremental  load  methods  described  in  [l]  and  [2]  should  be 
employed . 


VIII.  CONCLUSIONS 


This  report  has  dealt  primarily  with  the  additional 
mathematics  required  to  incorporate  elasticity  in  the 
deviatoric  stresses  developed  in  flowing  polymer  melts. 
Implementation  of  the  equations  within  an  existing  Finite 
Element  Computer  routine  was  then  shown.  From  these  analy¬ 
ses  we  can  make  the  following  conclusions: 

•  The  direct  Picard  Iteration  Converges  within  a 
radius  of  Ws£0.01. 

•  For  cross  channel  flow  and  entry  flow  "creeping" 
solutions  are  very  accurate  for  typical  polymer  extrusion 
Reynold's  numbers  (Re<0.4). 

•  For  the  Weissenberg  numbers  which  yielded  conver¬ 
gence,  no  appreciable  effects  on  the  flow  were  noted. 

•  The  programming  technique  of  passing  data  between 
elements  by  common  memory  appeared  to  be  effective. 

•  When  convergence  was  achieved,  the  calculated  values 
of  elastic  stresses  were  consistent  and  reasonable. 

•  The  penalty  method  of  incompressible  flow  appears 
to  yield  good  results  for  viscoelastic  fluids. 

•  The  radius  of  convergence  was  consistent  with  previous 
finite  element  calculations. 

•  The  radius  of  convergence  can  certainly  be  improved 
by  finite  difference  calculations  as  evidenced  by  Perera  [27]. 

•  Without  improvement,  the  only  computer  options  which 
should  be  used  in  evaluating  polymer  fluids  are  Newtonian  and 
power  law  viscous  (isothermal  and  adiabatic) . 


•  Techniques  o£  improving  the  viscoelastic  model  have 
been  proposed  which  offer  great  potential. 

e  For  24  element  problems,  the  computer  cost  for  runs 
requiring  30  iterations  was  $100.00. 
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IX.  RECOMMENDATIONS 


It  is  felt  that  the  work  performed  in  this  study  offers 
potential  for  useful  follow-on  effort.  In  particular,  there 
are  three  areas  of  development.  First,  the  analysis  of  the 
complete  flow  problem  is  vital.  While  the  gyroscope  fabrica¬ 
tion  is  new,  the  need  for  numerical  evaluation  in  the  molding 
process  is  not.  The  work  at  Cornell  [11]  demonstrates  this 
fact.  In  that  effort,  the  various  regions  of  flow  are  being 
tied  together.  A  similar  approach  is  required  for  the  finite 
element  modeling.  A  model  which  connects  the  flow  within 
and  out  of  the  extruder,  through  the  various  conduits,  and 
into  the  mold  cavity  is  an  important  development  which  should 
be  pursued. 

Direct  extensions  of  the  work  addressed  in  this  study 
are  also  important.  The  approach  should  be:  (i)  ignore 
stress  gradients  in  the  constitutive  equation  and  conduct 
direct  calculations,  (ii)  add  stress  gradients  along  with 
continuation  solution  methods  of  non-linear  equations.  Even 
if  future  work  with  constitutive  models  which  include  stress 
gradients  are  unsuccessful,  it  is  felt  that  the  equation 
with  some  elastic  stresses  will  be  a  big  improvement  over 
Newtonian  or  power  law  fluids. 
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Finally  as  efforts  one  and  two  above  progress ,  there  is 
a  need  to  conduct  rheology  experiments  which  will  determine 
properties  of  the  fiber-filled  polymers  being  used  in  the 
gyroscope  fabrication.  These  data  are  required  to  correlate 
with  the  velocities  and  stresses  predicted  by  finite  element 
equations . 

The  three  categories  are  listed  below: 

•  Model  complete  flow  history  from  extruder  to 
mold  cavity. 

e  Refine  viscoelastic  model. 

•  Conduct  rheology  experiments  of  appropriate  polymeric 


materials. 
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Figure  1  Piece  Part  Assembly  of  a  Plastic  Gyroscope 
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Figure  4  Non-Uniform  Molecule  Mesh  for  Solving  Stress 
Gradients 
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o  NODES 
®  GAUSS  POINTS 

EQUATION  IV-18  USED  AT  ®i<j 

EQUATION  IV-23  USED  TO  CALCULATE  0gN  P  AT  © 

EQUATION  IV-19  MODIFIED  AT  ®i'i+1  AS  FOLLOWS: 

xi,J+1  =  x©  yij+1  =  y©  aU+1  =  o® 


(Note:  Superscripts  are  indexed  at  each  Gauss  point  so  that 
whereas  it  is  x''i+1  referred  to  Gauss  point  4) 


x©isxU+2 


referred  to  Gauss  point  1 


Figure  5  Calculation  of  Elastic  Stresses  at  Gauss  Points 
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Figure  6  Iteration  Schemes  for  the  Solution  of  Creeping, 

Viscoelastic  Finite  Element  Equations:  (a)  Nested 
Iteration  (b)  Combined  Iteration 
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Figure  8  Flow  Geometries  and  Boundary  Conditions: 

(a)  Cross  Channel  Flow  (b)  Plane  Couette  Flow 
(c)  Entry  Flow  (d)  Step  Flow 


Velocity  Flow  Field  for  Linear  Cross  Channel 
Flow  (Vectors  scaled  relative  to  Un=100cm/sec) 


Vg  *  100  cm/s 


Figure  10  Velocity  Comparisons  of  9  and  8  Node  Elements: 

(a)  Element  15  (b)  Element  9  (c)  Element  7 
(Dashed  Arrows  are  8  Node  Elements) 
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(c) 

Figure  12  Fully  Developed  Flow  Behavior  of  Viscoelastic 

Fluid  Entering  and  Leaving  a  Contracting  Channel: 

(a)  Vortex  angle  3  (after  White  [33])  (b)  3  vs  Ws 

(after  White  [33])  (c)  Finite  Difference  Calculation 
for  Ws=0.6  (after  Perera  [27]) 
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Figure  15  Velocities  and  Pressures  for  Newtonian  Entry 
Flow  (Curve  for  velocities  is  schematic  only) 


study) 


Figure  17  Solutions  to  Nonlinear  Equations  (a)  Picard 
Iteration  (b)  Newton-Raphson  Iteration 


RUN  NO. 


GEOMETRY 


TYPE 


CONVERGENCE 
(YES  OR  NO) 


COST  (S) 


1 

CROSS CHANNEL 

LINEAR 

18-9  NODE  ELEM. 

- 

2.00 

2 

CROSS  CHANNEL 

LINEAR 

18-8  NODE  ELEM. 

- 

2.00 

3 

CROSS  CHANNEL 

LINEAR 

72-8  NODE  ELEM. 

- 

3.00 

4 

CROSS CHANNEL 

CONVECTION  (Re  =  0.4) 

18-9  NODE  ELEM. 

YES 

3.50 

5 

CROSS CHANNEL 

VISCOELASTIC  <WS=0.1) 

18-9  NODE  ELEM. 

NO 

10.00 

6 

CROSS  CHANNEL 

VISCOELASTIC  (WS  -  0.02) 
18-9  NODE  ELEM. 

YES 

21.35 

7 

CROSS  CHANNEL 

VISCOELASTIC  (WS  =  0.06) 
18-9  NODE  ELEM. 

NO 

32.96 

8 

PLANE  COUETTE 

LINEAR 

18-9  NODE  ELEM. 

- 

2.00 

9 

PLANE  COUETTE 

VISCOELASTIC  (WS  =  0.06) 
18-9  NODE  ELEM. 

NO 

12.47 

10 

PLANE  COUETTE 

VISCOELASTIC  (WS=0.02) 
18-9  NODE  ELEM. 

YES 

23.41 

11 

ENTRY 

LINEAR 

24-9  NODE  ELEM. 

— 

2.00 

12 

ENTRY 

VISCOELASTIC  (WS  =  0.01) 
24-9  NODE  ELEM. 

TENDING  AT 

30  ITERATIONS 

87.77 

13 

ENTRY 

VISCOELASTIC  (WS  =  0.001 ) 
24-9  NODE  ELEM. 

YES 

77.00 

14 

ENTRY 

VISCOELASTIC  (WS=  0.03) 
24-9  NODE  ELEM. 

NO 

37.39 

15 

STEP 

LINEAR 

30-9  NODE  ELEM. 

- 

2.50 

16 

STEP 

CONVECTION  (Re  =  0.4) 

30-9  NODE  ELEM. 

YES 

21.46 

17 

STEP 

VISCOELASTIC  (WS-0.01) 
30-9  NODE  ELEM. 

YES 

115.00 

18 

STEP 

VISCOELASTIC  (WS  =  0.001 ) 
30-0  NODE  ELEM. 

YES 

79.12 

19 

STEP 

VISCOELASTIC  (WS-0.03) 
30-9  NODE  ELEM. 

NO 

50.71 

Table  Computer  Run  Matrix 
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APPENDIX  1 


Derivation  of  Elastic  Stress  Gradient  Expressions 


From  figure  3  we  can  write  the  Taylor  series  approximations 
for  7 o  as: 


Forward  Difference:  a*"*’3''3  =  a3-' 3  +  |^ii,jAxf  +  jAyf  + 


i  — li  ^x2  +  T  — li  ,Ay|  +  ... 

7  3x2  X'3  f  23y2'i,D  f 


°i>j+i  -  »i>j  ♦  Mil.:**; ♦  + 

_*2,  1  3^a  : 

3y 


1  i!£|  A  *f+  1  2i£|  ,y*2+  ... 

2  .-z!i,3  f  2  .,2*1,3 


3x 


Backward  Difference:  1,3  =  a3"'3  -  Ix^i  jAxb  “  fy^i,jAyb  + 


+  •  •  • 


1  i!®|  Ax2  +  t  — L  ,Ay2 
7  3x2  1,3  b  1  3y2  1,3  b 

i!£|  iAy*2+  ... 

3x2  1,3  ”  1  3y2  1,3  b 


where : 


Axf  =  xi+1'3  -  x1'3  ,  Ayf  =  yi+1'3  -  yi/3  , 
Ax!  -  x1'^1  -  x1'3  ,  Ay!  =  yi'3+1  -  y*'3  , 


Axb  =  xi,j  -  xi_1,3  ,  Ayh  =  y1'3  -  y1"1'3  , 


:f 

'b 


83 


and 


s*b  -  1‘i'3  -  xi'j'1  -  *yb  ■  yi,)  -  y1,1’1 


Subtracting  the  first  and  second  equations  of  the 
backward  differences  from  the  respective  forward  differences: 


,i+1 • 3  _  ai“lfj  =  3£ 


+  lxb>  +  lfii,j(4yf +  ^b1  + 


1  32a  I 


3x 


Ay2)  +  ...  0 (A3) 


°1#3+1  -  =  !£li,j(Axf  +  Axb?  +  If 'i,  j  (A^f  +  A^>  + 


1  32a i  *2  .  *2  .  1  32o i  .. 

2  — -L  ^(Axf  -  Ax,  )  +2  — .  (Ayf 

i  3x2  f  2  3y2  1,3  f 


*2 

-  Ayb  )  +  ...  0 ( A  ) 


Assuming  that  all  differences  of  the  intervals  squared 
are  infinitesimal  (zero  for  the  uniform  mesh  case)  and 
solving  for  the  gradients  we  have  in  matrix  form: 


(Axf  +  Axb)  (Ayf  +  Ayb) 
(Axj  +  Ax£)  (Ay*  +  Ay*) 


r  3a  1 
TZ 

’  ai+1,j  -  o1-1^ 

3a 

_i» j+1 

_  i,  j-1 

UyJ 

U 

-  a  J 

We  can  use  Cramer's  rale  for  the  solution  since  the 
determinant  of  the  coefficients  of  the  gradients  can  never 
vanish.  Therefore: 
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D  X  Ol>i 

rt)  CO  CO  |*T> 


(a1+1,j  -  a1“1,j)(Ay*  +  Ay*)  -  (a1,j+1  -  a1'j_1)(Ayf  +  Ayb) 
(Axf  +  Axfc)  ( Ay^  +  Ayb)  -  (Axf  +  Ax£)  (Ayf  +  Ayb) 

(o1/j+1  -  ai'j"1)(Axf  +  Axb)  -  (oi+1'j  -  a1_1'j)(Ax*  +  Ax*) 

*  *  It  it 

(Axf  +  Axb) (Ayf  +  Ayb)  -  (Axf  +  Ax^ (Ayf  +  Ayfa) 


When  substitutions  are  made  for  the  A  terms  we  obtain 
equations  IV. 19. 
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APPENDIX  2 


Calculation  of  the  Global  Second  Derivatives 

A  subroutine  ESHAP  was  written  to  calculate  the  global 
second  derivatives  of  the  velocity  vector.  For  a  nine-node 
Lagragian  isoparametric  element  the  trial  functions  are: 


Ni 

-  r)  (s2 

-  s) 

n2 

+  r) (s2 

-  s) 

n3 

■i«** 

+  r) (s2 

+  s) 

N, 

- 

-  r) (s2 

+  s) 

Ns 

-  1) (s2 

-  s) 

Ns 

+  r)  (s2 

-  1) 

n7 

II 

1 

to|  1— • 
h" 

-  1) (s2 

+  s) 

Nb 

-  ^  (r2 
~  1{  r 

-  r) (s2 

-  1) 

N, 

=  (r2 

-1)  (s2  - 

1) 

We  can  form  the  following  table: 
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AD-A106  740  AIR  FORCE  INST  OF  TECH  KRX6HT-PATTERS0N  AFB  OH  F/6  20/11 

A  FINITE  ELEMENT  MODEL  OF  A  WHITE-METZNER  VISCOELASTIC  POLYMER  — ETC(U1 
FE8  81  B  R  COLLINS 

UNCLASSIFIED  AFIT-C I-A1-31T _  HL 


Where  the  Jacobian 


J  = 


3x  dx 
5r  3s 


has  been  used.  All  the  terms  in  this  equation  are  available 
at  the  Gauss  points  e.g. 


G.P . 


Xi 


G.P. 


where  are  the  x  coordinates  of  node  i. 

We  can  then  solve  for  the  global  second  derivatives  according 


APPENDIX  3 


Listing  of  New  Subroutines 


1 .  ELMT05 

2.  ELMT06 

3.  ESHAP 

4 .  PFORM 

5.  CMATRX 

6 .  FPSIG 


t 
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1 


BRC4066  (FOREGROUND)2  OUTPUT  FROM  TSO  XPRINT 
AT  10-* 92 : 40  ON  H/07/60  -  BRC4066.ELMT0S.FORT 

SUBROUTINE  ELMT05(  0  .  UL  .  XL  .  IX  .  TL  .  S  .  P  . NOF ,NOH . NST , ISM ) ° 

ELKT05  — — . ******!!S2^SS«2 

C*»***»**«**«»******<m»*  »**»ihhmmmmmmmhhmmnhhm«**»mmmh***OMCO0§O 

C  A  GENERAL  PENALTY  ELEMENT  FOB  INCOMPRESSIBLE  FLUID  FLOM  SJmOOBO 

#S85 

1  XM( 4 ,4 1/2*1 .DO ,4*0 .00 • 2*1 .DO #4*0.00 , 5*1 .00 . 3*0 .DO , 3*1 .00,3*0.00/  00000110 

JmMON*/raATA^oiHEAofEOI.NUMNP,NUMEL,«»t1AT,NEN.NE(>.IPR  222oOx’o 

COMMON  /ELDATA/DM.N.MA.MOT.IEl.NEL  00000150 

COMMON  /TAYLR/ESIG1(4,2,50)»ESIG2(4»2,50),ESIG3(4»2»50)»YY(4#2»5G  00000160 


1  ),EUSl(4,2,50),ELAS2(*.2.50),ELAS3(4.2.SO),EeSIG(4,3,50) 
DIMENSION  0< 30 )#UL(NDF ,ll#XL(NOM.l )#IXl 1 1 >TLI 1), 

1  S(NST,1).P{1).3HP(3.9I.SG<4).TG(9),KG(9),  .... 

2  SIG(7).EPS(6).BSIG(3).XX(3).BI10).DB(6,3).BT0BI3,3># 

3  BU(6),XMT3(3I.XMTST(3I.PEN(3,3),0U(S).DLTEE(3). 

4  V(2)  .DV(2.2I  ,XN(2.2)  ,A0VEC(2»2)  ,CADVEC(2»2) 

5  #  A1TERC  3,3), ESHP( 3 • 9 ) #00V( 3,2) 

DATA  PI/3. 141592653600/ 

IF  (1SU.EQ.1)  GO  TO  1 
ITYPE  =  0(30) 

L  *  0(2«) 

RHO  *  0(27) 

XLAM  s  0(26) 

XMU  =  0(25) 

XK  «  0(24) 

C  *  0(23) 

N1  «  0(20) 

HEAT  r  0119) 

LLB  3  LB( ITYPE) 

K2  3  0(18) 

N3  3  0(17) 

8  3  0(16) 

P4  3  0(15) 


C  BRANCH  TO  CORRECT  ARRAY  PROCESSOR 
C 

^  GO  TO  (1.2. 3, 3. 5,3,3). ISN 


ISM  «  1!  READ  MATERIAL  PROPERTIES,  DEVELOP 
DIAGONAL-STORAGE  0  MATRIX 


1  CALL  OFMTRX(D) 
LINT  »  0 

RETURN 


00000165 

00000170 

00000180 

00000190 

60000200 

00000210 

00000215 

00000220 

00000210 

00000240 

00000250 

00000260 

00000270 

00000280 

00000290 

00000300 

00000310 

00000320 

00000330 

00000340 

00000350 

00000360 

00000370 

00000380 

00000390 

00000400 

00000410 

00000420 

00000430 

00000440 

00000450 

00000460 

00000470 

00000480 

00000490 

00000500 

00000510 

00000520 

00000530 

00000540 
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II 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

c 

c 


301 

302 

c 

c 

c 


31 


32 

C 

c 

c 


320 

325 


9 

-C 

c 


2  RETURN 


ISM  *  j:  FORM  ELEMENT  STIFFNESS  MATRIX 


3  CONTINUE 

LOOP  OVER  GAUSS  INTEGRATION  POINTS 
COMPUTE  UNSTMMETRIC  STIFFNESS  MATRIX 

IF  I L**NDM  .HE.  LINT)  CALL  PGAUS3  ( L. LINT.SG.TG.MG ) 

DO  33  LL*1,LINT 

CALL  SHAPE  t SGI LL). TGI LL).XL.SHP,X3J.NDH,NEL. IX, .FALSE. ) 
WGT*X3J«MS( LL) 

COMPUTE  RADIUS  FOR  AXISVMMETRIC  CASE 

IF  (ITYPE.NE.3)  GO  TO  302 
RR=0.00 

00  301  1*1. NEL 

PR=RR*SHP<  3. 1  )»XU  1.1) 

CONTINUE 

MGT=MGT*2.D0»PI*RR 

CONTINUE 

COMPUTE  COORDINATES.  VELOCITIES  AND  GRADIENTS  FOR  CONVECTIVE 

DO  32  1*1 .NON 
XXI 11*0.00 
V(I)*0.00 
DO  31  K*1.NEL 

XX(  I  )=XX(  I  )+SHP(  3.K  >»XU  I  .K ) 

V(  I  )=V(  I  )+SHP<  3  .K  )»UL(  I  ,K ) 

CONTINUE 

YY(LL.I.N)  *  XXII) 

DO  32  J*1.N0M 
OV(I,J)=0.00 
DO  32  K=1,NEL 

DV(I,J)=DV(I,J)+SHP(J,K)»UL<I,K) 

COMPUTE  NONLINEAR  VISCOSITY  CORRECTION 
XNLNP*1.00 

IF  (P4.EQ.1.)  GO  TO  325 

Al*2 . 00*( 0V( 1 , 1 )«*2«0V< 2 , 2 }*#2 )♦ ( 0V( 1 , 2 )+DV( 2 , 1 ) )**2 
IF  C ITYPE.NE.3)  GO  TO  320 
A1*A1+2.00*(V( 1 )/XX( 1 ) )**2 
XNlMR=XNLNS/< l.DO+Al**< ( 1 .D0-P4 )/2 .00 ) ) 

VISUM  *  0.00 

IF  (G.EQ.O.DO )  GO  TO  9 

VISUM  *  XNLHR«XMU/G+VISLAM 

IF  (ISM. EG. 6. OR. ISM. EG. 4. OR. ISM. EQ. 7)  GO  TO  47 

LOOP  OVER  COLUMNS,  FORMING  OB,  MT*B.  AND  <DEL.(NU)T)T»N 

DO  46  J*1,NEL 

CALL  BMATRX(B,J,ITYPE,SHP,RR) 


00000550 
00000560 
CC0IC570 
00000500 
00000590 
00000600 
00300610 
00000620 
00000630 
00000640 
00000660 
000C0670 
00000600 
00000601 
00000690 
00000700 
00000710 
00000720 
00000730 
00000740 
00000750 
00000760 
00000770 
00003760 
00000790 
00000600 
00000610 
00000820 
00000630 
TERM  00000640 
00000S50 
00000860 
00000870 
00000880 
00000890 
00000900 
00300910 
00000920 
00000925 
00000930 
00000940 
00000950 
00000960 
00000970 
00000980 
00000990 
00001000 
00001010 
00301020 
00001030 
00001040 
00001050 
00001070 
00001072 
00001076 
00001090 
00001100 
00001110 
00001120 
00001130 


I 
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57 


C 

C 

c 

c 

I 


c 

c 

c 


c 

c 

c 


45 

46 

47 

C 

C 

C 

C 

C 


C 

C 

C 
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CALL  VMU10F(D,LLB.B.N0M,LLB.DB.6I 

CALL  VMULFF(XMT<ITYPE.1),8,1,LLB,NDM,4,LIB,XMTB,1,IER( 

DO  37  IDEXS1 >KCM 
00  37  JDEX-1 ,NDM 

IF  IIDEX.EQ.JOEX)  XNf IDEX, JDEX )sSHP( 3. J ) 

IF  (IDEX. NE. JDEX)  XN< IDEX, JOEX )=0 .DO 

CALL  VMULFF<DV,XN,NDH,NDM,NDM,NOM,NOM.AOVEC,NDM,IER> 

JJ*( J-l )*NDF+1 

LOOP  OVER  ROUS.  FORMING  BT«(DB).(MTB)T«MTB.  AND  NT(DEL.(NU)T)T*N 
DO  45  1*1. NEL 

CALL  BMATRXfB.X.XTYPE.SHP.RR ) 

CALL  VMULFMf B.OB.LLB.NOM.NDM, LLB.6.BT0B , 3. IER ) 

CALL  VMULFPfXMTf ITYPE.I ) >B. 1. LLB.N0M.4. LLB.XMTBT.l.IER ) 

CALL  '  VWLFMf XMTBT  .XMTB .  1 ,  NOM .  NDM ,  1 , 1  >  PEN ,  3 . 1 ER ) 

CALL  VMULFMf XN.AOVEC .NDM. NDM, NOM. NDM. NOH. CADVEC .NON, IER) 

II*( I-l )*NDF*1 

ADD  TO  ELEMENT  STIFFNESS  MATRIX  SINST.NST) 

CALL  MXAD0iS(II,JJ),NST,BTDB,3,N0M,N0M,WGT«XNLNR> 

CALL  MX  ADD  ( S(  II ,  J  J 1  .NST ,  FEN ,  3  .NDM  ,NDM  ,t:GT*XLAM ) 

CALL  MXADD  ( Sf II , J J ) ,NST , CAO VEC , NDM , NDM , NDM ,  L'GT*RHO ) 

ADD  THERMAL  STIFFNESS 

IF  (Nl.EQ.l )  A2=XK*00T(SHP( 1 ,1 ) ,SHP( 1 , J ) ,NDM ) 

IF  (Nl.EQ.l)  S(II+N0H,JJ*N0M)=SfII*NDM,JJ+N0M)+A2*WGT 

CONTINUE 

CONTINUE 

IF  ( I5U.EQ.3)  GO  TO  65 
CONTINUE 

IF  ( ISU. EQ.4.GR . ISW.EQ.6 )  GO  TO  60 

CALCULATE  ESIG( LL.2.N) :  ELASTIC  STRESS  AT  K  ♦  1  ITERATION 
SET  UP  A  MATRIX  FOR  PLANE  FLOW 


AITER(l.l) 
AITERf  2.1 ) 
AITER( 3,1 ) 
AITER( 1,2  ) 
AITERt  2,2 ) 
AITERf  3.2 ) 
AITERf 1,3) 
AITERf 2,3) 
AITERf 3,3) 


=  DV( 1,1 )*2.D0 
=  O.DO 
=  DVf 2,1 ) 

=  O.DO 

=  DVf  2  >2 )*2 .DO 

*  0V(1,2) 

*  DVf 1,2 >*2.00 

*  DVf 2 ,1  )*2.DQ 

=  DVf  1,1)  ♦  DVf 2,2) 


COMPUTE  VISCOUS  STRESSES  AT  GaUSS  POINTS:  SIG  =  0*BU 


SIG(l)  *  XMU*2 .DO*OV( 1,1 )*XNLNR 
SIGf 2 )  *  XMU»2.D0*DV(2,2)*XNLNR 
3IG(3)  *  XMU*(DV( 1,2 )  ♦  DVf 2,1 ) )«XNLNR 
SIGf 7)  =  XLAM*( DVf 1 ,1 )  ♦  DV(2,2)) 

IFf ITYPE.NE . 3 )  GO  TO  61 
SIGf 4)  *  SIGf 3) 

SIGf 3)  *  ( V( 1  )/XX( 1 ) )*XMU*XNLNR*2.D0 
SIGf 7)  =  SIGf 7)  4  XLAM*( V( 1 )/XX( 1 ) ) 


00001140 

0C0C1150 

00001160 

00001170 

00001100 

00001190 

00001200 

00001210 

00001220 

00001230 

00001240 

00001250 

00001260 

00001270 

00001280 

00001290 

00001300 

00001310 

00001320 

00001330 

00001340 

00001350 

00001360 

00001370 

00001380 

00001390 

00001400 

00001410 

00001420 

00001430 

00001470 

00001480 

000014S5 

00001487 

00001497 

00001507 

00001517 

00001527 

00001537 

00001547 

00001557 

00001567 

00001577 

00001587 

00001597 

00001607 

00001617 

00001627 

00001637 

00001647 

00001657 

00001667 

00001677 

00001667 

00001697 

00001707 

00001717 

00001727 

00081737 

00001742 


II 


94 


61 

C 

C 

C 

c 

c 

c 


21 


22 

C 

c 

c 


c 

c 

c 


c 


IFt ISW. EQ.4.CR. ISW.EQ.6 )  GO  TO  62 

CALCULATE  VISCOUS  STRESS  GRADIENT  (DEL(SIGMA) ) 

CALL  ESHAP(SG(  LD.TGt  LL1  .XL.ESHP.NDM.NEL.IX) 

FCSM  CONVECTION  DERIVATIVE  OF  STRESS)  ONLY  2D  FLOW 


DO  21  1=1.2 
DO  21  J-1,3 
OOV(J.I)  =  0.00 
DO  22  K=1.NEL 
ODV(l.l)  s  ODV(l.l) 
D0V(1,2>  =  ODV( 1.2 ) 
DOVt  2 . 1 )  =  D0V(2.1) 
DDVI2.2)  =  DOV( 2.21 
DDVC3.1)  =  0DV(3.1) 
1  l.K)> 

DDV( 3.2 )  =  ODVt  3.2 ) 
1  l.K) ) 

CONTINUE 


♦  2.D0»XMU*XNLNR»ESHPtl,K)*ULtl,K) 

♦  2 .DO»XHU*XNLNR*ESHP< 3.K )*UL< l.K) 

♦  2  .D0*XMUl<XNLN3*ESHPt 3 .K )»Ul( 2 ,K ) 

+  2 .DO*XMU*XNLNR*ESHP<  2  ,K  )»UU  2  ,K ) 

♦  XHU*XNLN9*<  ESHPt  1  ,K  >«UL(  2  ,K  UESHPt  3.K  >»ULt 

♦  XMU»XNLNR»( ESHPt  3 ,K XULt 2 ,K KESHPt 2 ,K >*ULt 


SOLVE  ESIGt  LL.2 ,N) t  ONLY  2D  FLOU 


ESIGlt  LL.2 ,N)  =  VISLAM»l ( AITERt 1.1 )*t SIGt 1 )  +  ESIGlt LL.l.N) )4 

1  AITERt 1 ,2 )*( SIGt  2  HESIG2I LL. 1 .N) )♦ AITERt 1 ,3  >*l  SIGt  3  KESIG3I 

2  LL.l,N)))-tV(l )*(COV( 1.1 l+ELASlf LL.l >N)  )+Vt  2 )«t  DOVt 1 ,2  HELASlf 

3  LL.2.N) ) ) ) 

E3ZC2( LL.2.N)  s  VISLAM«tf AITERt 2,1 )»I SIGt 1 HESIGlt LL.l, N ) )♦ 

1  AITERt  2 ,2  )*t  SIGt  2  )+ESIG2<  LL .  1  .N )  )♦ AITERt  2 ,3  >*l  SIGt  3 HESIG3C 

2  LL,l.N)))-tVtl )*t DOVt  2 .1  HELAS2I LL.l >N) )+V( 2 )*l DDVl 2.2 1+ELAS2I 

3  LL.2.N) ) ) ) 

ESIG3t  LL.2 ,N)  =  VISLAM»< ( AITERt  3,1 )*( SIGt 1 1+ESIGlt  LL.l.N))* 

1  AITERt  3,2 )*t SIGt  2 )+ESIG2f LL.l ,N) )+AITERt  3.3 )»f  SIGt  3 )+ESIG3t 

2  LL.l >N ) ) )-( Vt 1  )»< ODVt  3.1 )+ELAS3f LL.l ,N ) )*V( 2  )#( DDVt  3,2 )*ELA53t 

3  LL.2.N)))) 


UPDATE  BOUNDARY  STRESSES  BOSIGtNOOE, DIRECTION, EU1T.  NO.) 


IF  IG.EQ.O.OO)  GO  TO  65 

BOSIG(LL.l.N)  S  ESIGlt  LL.2.N)  ♦  ELASK  LL,l,N)*(XLf  1,LL>- 
1  YYt  LL.l.N) )  ♦  ELASK  LL.2 ,N )*(  XLt  2  ,  LL  )-YY(  LL.2  ,N ) ) 

BOSIGt  LL.2.N)  =  ESIG2t  LL.2.N)  ♦  ELAS21  LL.l ,N>*( XLt  1, LL>- 
1  YYt  LL.l.N) )  ♦  ELAS2( LL.2 ,N )»( XLt  2 , LL)-YYt  LL.2 ,N) ) 

BOSIGt  LL.3.N)  =  ESIG3( LL.2.N)  ♦  ELAS3t LL.l ,N )*t XLt 1 , LL >- 
1  YYt  LL.l.N) )  ♦  ELAS3( LL.2  >N)»t  XLt  2 , LL)>YY( LL.2.N) ) 

GO  TO  65 

PRINT  STRESSES  IF  ISW=4,  OTHERWISE  BRANCH  TO  COMPUTE 
UNBALANCED  FCSCE  VECTOR 

IF  IISW.EG.6)  GO  TO  66 

XMAX  *  DMAXltOABSf XLt 1.4 )-XLt 1.1)) .OABSt  XLt  1.3)-XLt  1,2 ) ) , 

1  OABSt  XLt  2 .4 )  —XLt  2,1))  .DABSt  XLt  2 ,3  )-XL(  2,2))) 

SIGt  5)  =  (RHO/t  XMU*XNLNR )  )*OSOST(  Vt  1  )**»2*Vt  2  )»*2  )»XMAX 

SIGt  6 )  =  .  ISLAM«DStJRTt  Vt  1  )»*2*V(  2  )**2  l/XMAX 

CALL  FPSIGt XX, ESIGlt  LL.2 ,N) >ESIG2( LL.2 .N ). ESIGlt LL.2 ,N) ,SIG, 

1  ITYPE.NDF) 

GO  TO  65 


H 


00001747 

00001757 

00001767 

00001777 

00001787 

00001797 

00001807 

00001817 

00001827 

00001837 

00001847 

00001857 

00001867 

00001877 

00001S87 

00001697 

00001907 

00001917 

00001927 

00C01937 

00001947 

00001957 

00001967 

00001977 

00001987 

00001997 

00002007 

00002017 

00002027 

00002037 

00002047 

00002057 

00002067 

00002077 

00002037 

00002097 


00002107 

00002117 

00002127 

00002137 

00002147 

00002157 

00002158 

00002159 

00002160 

00002161 

00002167 

00002177 

00002187 

00002197 
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LOOP  OVER  NODES  TO  COMPUTE  UNBALANCED  FORCE  VECTOR: 

P  =  PI  -  BT*5IG  -  NT»ELASCLL,2,M)-RHO*tNTCDEL.(KUmTN 


00002207 

03002217 

00002227 

00002237 

*0002247 


COMPUTE  UNBALANCED  TEMPERATURE  VECTOR  00002237 

IF  (Nl.NE.l)  GO  TO  76  *0002247 

O  =  HEAT«( SIG( 1 )*DV( 1 , 1 l+SIGC 2 )»OV< 2 • 2 l+SIGC 3 )*CDVC 1 ■ 2 )+DV(  2 .1 ) ) )  00002257 
IF  CITYPE.EQ.3)  0  =  Q  ♦  HEAT*(0V< 1 ,2 HDVC 2 ,1 )  )*( SIG< 4 )-SIGC 3 ) )  00002262 

00  78  J=l,2  00002267 

OLTEE(J)  =  0.00  00002277 

DO  78  1=1, NEL  00002287 

DLTEECJ)  =DELTEE(  J>  ♦  SHPC  J,I  )«UL<  3,1 )  00002297 

DO  77  1=1, NEL  00002307 

II  =  CI-1HNOFU  00002317 

00002318 

CONVECTION  TERM  SAME  FOR  2D  AND  AXISYMMETRIC  FLOW  00002319 

00002320 

P(II)  =  P( II )  -  RHO»SHP( 3,1 )*C  VC  1 )*0VC 1,1 > ♦ Vt 2 )»0VC 1,2) )»WGT  00002321 

PCII+1)  =  P( 11*1 )-RHO*SHPC  3,1 )*<  VC  1 )»OVC  2 , 1 )+V( 2 )*OVt 2,2) )*MGT  00002322 
IF  ( ITYPE.EO. 3)  GO  TO  79  00002324 

PCII)  =  P< II )-(SHP( 1 ,1  )<H SIG( 1 l+SIGC  7) )*SHPC  2,1 )*SIG(  3 ) )«WGT  00002327 
PCII+1)  =  PCII+D-C SHPl  2 . I )»( SIG( 2  I +SIGC  7 )  l+SHPC  1,1  )*SIG( 3 )  )*WGT  00002337 
GO  TO  80  00002342 

PC  II I  =  PC  II )-( SHPt 1 ,1 )*C SIG( 1 )+SIG( 7) l+SHPC 3, 1 )*SIG( 3)  00C02343 

1  +SHPC2,I)*SIGC4))«WGT  00C02344 

PCII+1)  =  PCII+1)-CSHPC2.I)»CSIGC2)+SIG(7))+SHPC1,I)»SIG(4))«W3T  00002345 
IF  CK2.ED.3.0R.K2.EQ.4)  PCII)  =  PC  II l-CSHPC 3,1 )»C ELAS1C LL.1.N1+  00002347 

1  ELAS3CLL,2,N)))*HST  0C002357 

IF  CK2.EQ.3.0R.K2.EQ.4)  PCII+1)  =  PC  11*1 »  -  CSHPC3,I)«C  00002367 

1  ELAS2CLL,2,N)+ELAS3(LL,1,N)))»WGT  00002377 

IF  CN1.EQ.1)  A1  =  Q*SHP(3,I>  00002387 

IF  CN1.EQ.U  A2  =  XK«OOTCSHPC  1 ,1 )  .OLTEE.NOM)  00002397 

IF  CN1.EQ.1)  PCII+NDM)  =  PCII+NDM)  ♦  A1*WGT  -  A2»WGT  00002407 

CONTINUE  00002417 

CONTINUE  00002427 

RETURN  00002442 

END  00002447 

SUBROUTINE  ELMT06C  0  ,  UL  ,  XL  ,  IX  ,  TL  ,  S  ,  P  ,NDF ,NDM,NST,ISW)00COC01O 


**»»**»** 

C*«*»»**#***«»*********»X)*» 

c^*«***»**«**************** 

c 


C0000020 

«»«»««*««»**K«******»<tK»*»***00000030 

**tt****»***«)Hn,**»iHHi**»)HHnn»00000040 

**«*»****»*«*##»******»**»#*»0000C050 

00000055 


AN  ELEMENT  FOR  INTERPOLATING  DISPLACEMENT,  TEMPERATURE,  AND  STRESSOOOOOOfcO 
FOR  VISCOELASTICITY:  20  FLOW,  OLDROYD  DERIVATIVE  00000070 

00000080 

IMPLICIT  REAL*8C  A-H ,0-Z)  00000090 

REAL#8  XMTC 4,6 )/8»l.OO,2#O.D0,2»l.O0 , 12*0 .DO/,  00000100 

1  XNC6,4)/2*1.00,4XO.OO,2«1.DO,4«O.DO,3*1.DO,3»O.DO,3*1.DO,3*'O.DO/  00000110 


INTEGER  LBC 41/3,3,4,6/ 

COMMON  /COAT A/O  .HEAD  CEO),  NUMNP ,  NUMEL ,  NUMMAT  ,NEN  ,NE<J ,  IPR 

COMMON  /E  LD ATA/DM , N , M A , MOT , I E  L , NE  L 

DIMENSION  DC  30 ) ,UL( NOF.l ) .XLCNDM, 1 ) > IXC 1 ) ,TL( 1 ) • 

1  SCNST.l ) ,P( 1 ) ,SHPt  3,9 ) ,SG( 9 ) ,TGC  9) ,UG( 9 ) , 

2  SIGC  7) ,EPS( 6 ) ,BSIG( 3) ,XX( 3 ) ,B( 18 ) ,DB( 6 ,3 ) ,BTDBC  3,3) , 

3  BUC6) .XMTBC  3 ) .XMTBTC  3 ) ,PEH( 3,3 ) ,OUC  3 ) ,DLTEE( 3) , 

4  VC  2 ),DVC  2,2 ) ,XNC  3,3) ,ADVEC( 2,2 ) .CAOVECC  2,2 ) ,DDVC 3,2 ) ,C( 3,2 ) , 

5  XNTNC  3,3 ) ,BTC  2 ,3 ) , ADSIGC  3,2 ) ,CNC  3,2 ) ,XNTDB( 3,2 ) .XNTCNC  3,2 ) , 

6  XNTBTC 2,3) ,  CA0SIGC3.2) 

IF  (ISU.EQ.l)  GO  TO  1 


00000120 

00000130 

00000140 

00000150 

00000160 

00000170 

00000180 

00000190 

00000200 

00000205 

00200210 

00000220 


nnnw  w  non  o  oono  non  o  o  oono  non 


ITYPE  =  0(30) 
l  =  DCS) 

RHO  =  0(27) 

XLAM  =  0(26) 

XMU  =  0(26) 

XK  *  0(24) 

C9=  0(23)  . 

N1  *  0(20) 

HEAT  =  0(19) 

LLB  =  LB( ITYPE) 

K2  =  0(18) 

N3  =  0(17) 

6  =  0(16) 

P4  =  0(15) 

BRANCH  TO  CORRECT  ARRAY  PROCESSOR 
GO  TO  (1,2, 3, 3, 5, 3), ISM 

ISM  =  l:  REA0  MATERIAL  PROPERTIES,  DEVELOP 
DIAGONAL-STORAGE  0  MATRIX 

1  CALL  DrMTRX(O) 

LINT  =  0 

RETURN 

2  RETURN 

ISM  =  3:  FORM  ELEMENT  STIFFNES  MATRIX 

3  CONTINUE 

LOOP  OVER  GAUSS  INTEGRATION  POINTS 
COMPUTE  UNSYMMETRIC  STIFFNESS  MATRIX 

IF  (  L***NQM . NE .  LINT )  CALL  PGAUSS  (  L,  LINT.SG.TG.MG) 

DO  33  LL  =  1 , LINT 

CALL  SHAPE  ( SG( LL) ,TG( LL ) ,XL,SHP,XSJ,NOM,NEL,IX, .FALSE. ) 
MGT  s  XSJ*WS<  LL) 

COMPUTE  COOROINTAES, VELOCITIES,  STRESSES,  AND  GRADIENTS 

DO  32  1=1. NOM 
XX(I)=0.00 
V(I)=0.00 
DO  31  K=1 ,NEL 

XX(I>  =  XX(I>  ♦  SHP(3,K)*XL(I»K) 

V(I)  =  V(I)  ♦  SHP(  3,K  )«UL(  I  ,K ) 

1  CONTINUE 

00  32  J=1,NDM 
0V(I,J)=0.00 
DO  32  K=I,NEL 

2  OVtI.J)  =  0V(I,J)  +  SHP(  J  ,K  )*UL(  I  ,K ) 

COMPUTE  NONLINEAR  VISCOSITY  CORRECTION 

XNLNR  =  1.00 
IF  (P4.EQ.1. )  GO  TO  325 


00000230 

00000240 

00000250 

00000260 

00000270 

00000200 

00000290 

000003C0 

00000310 

00000320 

00000330 

00000340 

00000350 

00000360 

00000370 

00000380 

00000390 

00000400 

00000410 

03000420 

00000430 

00000440 

00000450 

00000460 

00000470 

00000480 

00000490 

00000500 

00000510 

00000520 

00000530 

00000540 

00000550 

00000560 

00000570 

00000530 

00000590 

00000600 

00000610 

00000620 

00000630 

00000640 

00000650 

00000660 

00000670 

00000680 

00000690 

00000700 

00000710 

00000720 

00000730 

00000740 

00000750 

00000760 

00000770 

00000780 

00000790 

00000300 

00003810 

00000820 


A1  =  2 .00*1 DVI 1 >1  )4**2*0V(  2 ,2 )»»2  HIDVI  1,2  )*DV(  2,1)  )**2 
XNLNR  5  XNIN3/11.00*A1«»I<1.03-P4)/2.00)) 

325  VISUM  =  0.00 

IF  IG.EQ.0.D0)  GO  TO  9 
VISUM  =  XNLNPWXMO'G  ♦  VISUM 
9  00  320  1=1,3 

SIGH)  =  0  .DO 
DO  320  J=X,2 
OOV(I.J)  =  0.00 
00  320  KH.NEL 

SIG1 1 )  =  SIG(I)  ♦  SHP!3,K)»ULINDF-3aI,K) 

320  ODVtI.J)  =  ODVII.J)  ♦  SHP( J  ,K  )*UK  NDF-3+I.K  > 

IF  ( I5W.EQ.4.0R. ISW.EQ.6)  GO  TO  97 
C 

C  LOOP  OVER  COLUMNS  FORMIN  MT*B,<OEL.<NU)T)T*N,  BT, 

C  DELI N*SIGMA )#N»  N.  08,  AND  CN 

C 

00  96  J=I ,NEL 

CALL  BMATRXIB, J,ITYPE ,SHP,RR ) 

CALL  CM.\TRXIC,J,SIG,SHP) 

CALL  VMULFF I XMTI ITTPE , 1 1 ,B . 1 , LLB ,NOF- 3 ,4 , LLB , XMTB , 1 , IER ) 

DO  37  IDEX=1,2 
00  37  JDSX=1,2 

37  ADVECIIOEX.JOEX)  =  OVI IOEX, JDEX)«SHPI 3,J) 

DO  41  I0EX=1,3 

DO  41  JDEXn.3 

IF  I IOEX. EQ. JDEX )  XNI IOEX, JOEX)  =  SHPI3.J) 

41  IF  I IOEX. ME. JDEX 1  XNIIOEX.JOEX)  =  0.00 
BTI1.1)  =  SHPll.J) 

BTC2.1)  =  0.00 
BT( 1 , 2 )  =  O.DO 
BTI2.2)  =  SHPI2.J) 

BTI1.3)  =  SHPI2.J) 

8TI2.3)  =  SHPll.J) 

DO  39  I0EX=1,3 
DO  39  JDEX=1.2 

AOSIGI  IOEX ,  JOEX )  =  SHPI 3,  J  i»OOVUOEX,  JOEX) 

39  CNI IOEX, JDEX  J  s  SHPI 3. J )*C< IDEX, JDEX) 

CALL  VHULDFIO, LLB,B,KDM,LLB,DB,6 ) 

JJ  =  IJ-1)*NDF  +1 
C 

C  LOOP  OVER  ROWS,  FORMING  I MTB )T*MTB»  NT! DEL. I NU )T)T*N,NT*BT , 
C  NT! OE LI N*SIGMA )*N ,  NT*N,  NT*OB,  AND  NT*CN 

C 

00  4S  1=1, NEL 
C 

CALL  BMATRXIB, I, ITTPE, SHP.RR) 

CALL  VMULFF! XMTHTYPE,l),B,l,lLB,NDF-3, 4, LLB, XMTBT.l, IER) 
CALL  VMULFMI XMTBT ,XMTB,1,NOM»NOM,1,1 ,PEN,3,IER ) 

DO  38  I0EX=1,2 
DO  38  J0EXn,2 

38  CAOVECI IOEX, JDEX)  =  AOVECI IOEX, JDEX)»SHPI 3,1 ) 

00  40  IDEX=1,2 

DO  40  JDEX=1,3 

XNTCNI JDEX, IOEX)  =  CNI JDEX, IOEX )*SHP( 3,1 ) 

XNTOBI JOEX, IOEX)  =  DBI JDEX, IOEX >»SHPt 3,1 ) 

XNTBTI IOEX, JOEX)  =  BTI IOEX, JOEX )*SHP( 3,1 ) 

40  CAOSIGI JDEX, IOEX)  =  AOSIGI JDEX, I0EX)*SHPI3, I) 

DO  42  IDEX=1,3 

DO  42  JDEX=1,3 


00000830 

00C0C840 

00000050 

OCOOOCfcO 

00000870 

00000880 

00000390 

00000900 

00000910 

00000920 

00000930 

00000940 

00000950 

00000960 

00000970 

00000980 

00000990 

00001000 

00001010 

00001020 

00001030 

00001040 

0C091050 

00001060 

00001070 

00001030 

00001090 

30001100 

00001110 

00001120 

00001130 

00001140 

00001150 

00001160 

00001170 

00301180 

00001190 

00001200 

00001210 

00001220 

00001230 

00001240 

00001250 

00001260 

00001270 

00001280 

00001290 

00001300 

00001310 

00001320 

00001330 

00001340 

00001350 

00001360 

00001370 

00001380 

00001390 

00001400 

00001410 

00001420 


u 
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42 

C 

c 

c 


45 

46 

47 

C 

C 

c 

c 


XNTNdOEX.JOEX)  =  XN<  IDEX,  JOEX  l«SHP<  3.1 ) 

II  =  (I-l)*!OF  ♦  1 

*00  TO  ELEMENT  STIFFNESS  MATRIX  S(NST,NST> 

CALL  MXAOO 1 S( II , J J ) ,NST ■ PEN , 3 . NOM , NOM , UGT»XLAM ) 

CALL  MXAOO! S( II, JJ).NST,CADVEC,NDM. NOM, NDH.N3T«RH0) 

CALL  MXAOO! S< II » JJ*NDF-2 ) >NST . XKTBT.2 .2.3.UGT ) 

CALL  MXAOO <  S( II*NOF -2 , J J ) ,NST , CAOSIG , 3 , 3 , 2 . WST»VISLAM ) 
CALL  MXAOO  ( S(  IItNDF-2 ,  J  J 1  ,NST .  XNTCN .  3  >  3  >  2 .  -k'ST«VISLAM ) 
CALL  MXAOO! S( II*NDF-2 . J J 1 ,NST , XNTDB .3.3.2. -NGT«XNLNR ) 
CALL  MXAOO! S< II+NDF-2, JJ*NDF-2 ) .NST.XNTN. 3.3.3. UST ) 

ADO  THERMAL  STIFFNESS 

IF(Nl.EQ.l)  A2  *  XK*OOT!  SHP!  1 , 1 )  ,SHP(  1,  J )  .NON) 

IF(Nl.EQ.l)  SI II*NOM, JJ+NOM)  =  S(  II+NDM,  JJ*N0MNA2»WGT 

CONTINUE 

CONTINUE 

CONTINUE 

IF  ( ISU.EQ.3)  GO  TO  65 

PRINT  STRESSES  IF  ISH=4. OTHERWISE  BRANCH  TO  COMPUTE 
UN3ALANCE0  FORCE  VECTOR 

IF  IISU.EQ.6)  60  TO  66 

XMAX  =  0MAX1IDABSIXLI 1.4 ) —XL! 1,1)) .DABS! XL!  1»3)-XL(  1 ,2 ) ) . 
1  DABS!  XL!  2 ,4 1-XLI  2.11)  .DABS!  XL!  2 , 3 )  -XL!  2.2))) 

SIS!  5)  =  !RH0/!XMU*XNLNR))«0S2RTiVll)**2*V!2)*»2)«XMAX 
SIGI6)  =  VI3LAM«0SQRT(V(1)»«2*V(2)»*2)/XMAX 
CALL  FPSIGI XX. 0.00. 0.00,0. 00. SIG.ITTPE.NDF) 

GO  TO  65 


C 

C 

c 

c 

c 

c 

66 


78 

76 


00001430 

00C01440 

000C1450 

00001460 

00001470 

00001480 

00001490 

00001500 

00001510 

00001520 

00001530 

00001540 

00001550 

00001560 

00001570 

00001580 

00001590 

00001600 

00001610 

00001620 

00001630 

00001640 

00001650 

03001660 

00001670 

00001680 

00001690 

00001700 

00001710 

00001720 

00001730 

03001740 

00001750 

00001760 

00001770 


LOOP  OVER  NODES  TO  COMPUTE  UNBALANCED  FORCE  VECTORS 
PI  =  PI  -  NTBT»SIGMA  -  RHO*NT( DEL. ( NU IT IT«N 
P2  =  P2  -  NT! DEL! SIGMA ) )*V  -  NT •SIGMA  ♦  NT*0*L*V  ♦  NT*HT»V, I»SI6MA00001 780 

00001790 

COMPUTE  UNBALANCED  TEMPERATURE  VECTOR  00001800 

IF  !N1.NE.1>  GO  TO  76  00001810 

Q  =  HEAT*! SIG! 1 )»OV! 1 .1 ) *SIG( 2 )*DV! 2,2))  00001820 

DO  78  J=l  ,2  00001830 

DLTEE! J )  =  0.00  00001840 

00  78  1*1. NEL  00001850 

DLTEE!  J)  =  DLTEE!  J)  ♦  SHPU.I  )»UL! 3.1 )  00001860 

DO  77  1=1, NEL  00301870 

II  =  II-1*<*NDF  ♦  1  00001880 

P(II>  =  P!II)-!RHO»!V!l)»OV!l,l)*VIE)»DVIl,2))*IDOVI  1.1)  00001C90 

1  +0DV! 3.2  I ) )*SHP( 3,1 )*K3T  -  SHP! 1.1 )*SIGI 7 >*M6T  00001900 

P!  11*1 )  =  P( 11*1 )-( RHO*! V! 1 )iDV! 2 ,1 )*V( 2 )»DVI 2,2) )»(OOVI 2,2 )  00001910 

1  ♦OOV!3,l)))*SHP!3.I)»WGT  -  SHP!  2.1  >*SI6!  7 )»KST  00001920 

IF  (Nl.EQ.l )  A1  =  Q»S!IP(3,I)  00001930 

IF  (Nl.EQ.l)  A2  =  XK»DDT! SHP! 1,1), DLTEE, NOM)  00001940 

IF  (Nl.EQ.l)  P< II+NOM )  =  P(II*NDM)  ♦  A1*U3T  -  A2«KST  00001950 

P( II*NDF-2 )  =  P( II+N3F-2 )-( VISLAM*! DDV! 1,1 )*V( 1 •♦OOVI 1 ,2 )»V( 2 ) )  00001960 

1  ♦  SIG! 1 )  -  2.DO»XMU*XNLNR»OV(1,1>  -  2.C0*VISLAM»  00001970 

2  <SIG(l)*0V(l,l)*5IG(3>«DV(l,2)))*SHP(3,I)»t!GT  00001980 

P(II»NDF-1)  =  P(II*NDF-1)-(VISLAM*(DDV(2,1)»V(1)*DDV(2,2)*V(2()  00001990 

1  ♦  SIG! 2 )  -  2.D0*XMU«XNLN3»DV(2,2)  -  .D0«VISLAM»  00002000 

2  (SIG(2)«OV(2,2>*SIG!3)»DV!2,l)))»SHP(3,I)»WGT  00002010 

P( II+NDF )  =  P( II+NOF )  -  (VISLAM«(DDV(3,1)*V(1)*DDV(3,2)«V(2))  00002020 


99 


i 


1  *815(3)  -  XMU»XNLNR»I0V(1,2)*0V<2,1))  -VISLAM» 

2  ( SIS<  2  >»0V(  1 , 2 )*SX61 3 ) *0V( 1 ,1  )*SIG( 1  )»0V( 2 , 1 ) 

3  +SIG<  3)»0V(2,2)  I  >«SHP<  3,I)*UGT 

77  CONTINUE 

65  CONTINUE 

33  CONTINUE 

S  RETURN 

END 

SUBROUTINE  ESHAPI SS,TT,X,ESHP,NDM,NEL,IX) 

C 

ESHAP  ■**»***#«#*»»****»***#*»»«(**** 
»»»****«»#**«**»*»#****»»*«*»* 

C 

implicit  real«sia-h,o-z) 

C  SHAPE  FUNCION  ROUTINE  FOR  9  NODE  QUADRILATERALS  FOR  SECOND  DER. 

C 

DIMENSION  ESHPI 3,1 ),X(NOH,l ) ,SHP( 3,9) ,1X1 1 ) >BI6( 3.3) >XS( 2.2 ) ■ 

1  EBIG( 3.3) >EXS( 3.2 ) >SX( 2.2) .TEMPI  3) 

DATA  S/0.500/,T/1.00/',R/2.00/ 

C 

C  FORM  9-NODE  QUADRILATERAL  SHAPE  FUNCTIONS  FOR  SECOND  DERIVATIVE 
ESHP(l.l)  =  Sk|TT*«2-TT> 

ESHPI2.1)  =  S*ISS*»2-SS) 

ESHP13.1)  =  S»»2*IR»SS-T)«IR»TT-T) 

ESHP(  1,2)  =  ESHPU.l) 

ESHPI 2,2 )  =  S»<SS»»2*SS) 

ESHPI3.2)  =  S**2«(R»SS*T)«(R*TT-T) 

ESHPI 1,3)  =  S»ITT»»2*TT) 

ESHPI 2,3)  =  ESHPI 2,2) 

ESHPI 3,3)  =  S»*2*(R*SS*T)«(R»TT*T) 

ESHPI  1,4)  =  ESHPI  1,3) 

ESHPI 2, 4)  s  ESHPI 2,1) 

ESHPI 3,4 )  =  S»»2»(R«SS-T)*IR*TT*T) 

ESHPI 1,5)  «  -R*cSHP( 1,2 ) 

ESHPI 2, 5)  =  T-S5*»2 
ESHPI 3,5)  =  SS»IT-R»TT) 

ESHPI 1,6 1  =  T-TT»*2 
ESHPI 2, 6  I  =  -R«ESHP{ 2,2  F 
ESHPI 3,6 )  *  -TT«(R*SS*T) 

ESHPI 1,7)  =  -R*ESHP(1,4) 

ESHPI 2,7)  =  ESHPI 2,5) 

ESHPI 3,7)  =  -SS*IR»TT*T) 

ESHPI  1,8)  =  ESHPI  1,6 ) 

ESHPI 2,8)  =  -R*E5HP( 2,1 ) 

ESHPI 3, S )  =  TT»IT-R«SS1 
ESHPI 1,9)  s  -R4ESHPI 1,6) 

ESHPI 2, 9)  =  *R»ESHPI 2,5) 

ESHPI 3.9)  =  R«*2*SS»TT 
C 

C  CONSTRUCT  BIG  MATRIX  AND  ITS  INVERSE 

C 

CALL  SHAPE! S3,TT,X,SHP,XSJ ,NDM,NEL, IX,. TRUE. ) 

DO  130  1*1 ,NOM 
00  130  3*1,2 
XSII.J)  *  0.00 
■  DO  130  K*1,NEL 

130  XSII.J)  *  XSII.J)  *  X(I,K)*SHP( J,K) 

BIG! 1,1)  *  XSI 1,1 )»»2 
BIG! 2,1 )  *  XS(2,1)«»2 


00002030 

00002040 

00302050 

00002055 

00002060 

00002070 

000020S0 

00002090 

00000010 

00000015 

00000020 

00000030 

00000040 

00000050 

00000060 

00000070 

ooooooso 

00000090 

00000100 

00000110 

00000120 

00000130 

00000140 

00000150 

00000160 

00000170 

00000180 

00000190 

00000200 

00000210 

00000220 

00000230 

00003240 

0C000250 

00000260 

00000270 

00000280 

00000290 

00000300 

00000310 

00000320 

00000330 

00000340 

00000350 

00000360 

00000370 


00000380 

00000390 

00000400 

00000410 

00000420 

00000430 

00000440 

000^0450 

00000460 

00000470 


1 


100 


1 


8I6(3>1  >  =  XS(1,1)*XS<2,1>  08000430 

BIG( 1>2)  =  XSI 1(2 )**2  00000490 

816(2.2)  =  XS(2, 21**2  OOOOC5CO 

816(3,2)  =  XS( 1.2 )*XS( 2,2 )  00000510 

816(1.3)  =  2 .D0*XS( 1,1 )*XSl 1,2)  00003520 

616(2,3)  =  2.00*XS( 2 , 1 )*XS( 2,2 )  00000530 

816(3,3)  =  XS(  1,1  )**XS(  2,2 )  ♦  XStl,2)*XSt 2,1)  00000540 

CALCULATE  DETERMINANT  OF  BIG  00000550 

DET  s  BIG( 1,1  )*( BIS( 2.2 )*BI6( 3,3) -BI6( 3,2 )*BIG( 2,3) )-BI6( 2,1)*  00000560 

1  (816(1,2 )«BXG( 3>3)-BIGll«3)*BIG(  3,2  >  H8IGI 3,1 >*(8X6(1 >2 >*8161  2 , J  >00000570 


2  -BIG( 2 , 2 )*B16( 1,3)) 

FORM  INVERSE 

EBIG(l.l)  =  ( BIG( 2,2 )*BIG( 3,3)-BIG( 3,2 )*BIG(  2,3) )/OET 
EBIG( 2,1)  *  -(BIG( 1 ,2 )»EIG( 3,3 )-BIG( 3,2 )‘D1G( 1 >3) )/OET 
EBIG( 3,1 )  s  (BIG( 1,2  )*BIG( 2 > 3 )-BIG( 2.2  )*6IG( 1,3) )/DET 
EBIG(1,2)  3  -(BIG( 2,1)*QIG( 3,3)-BIG( 3,1 )«015I 2,3 ) VCET 
E6IG( 2,2 )  3  ( BIG( 1,1 )*BI6( 3,3 )-BlG( 3,1 )*BIG( 2,3) )/DET 
EBIGI  3,2)  =  -I BIG(  1  ■  1  )*BIG(  2 ,3)-BI6l  2,1  )>tOIG(  1,3)  )/DET 
E8I6( 1,3)  3  ( EXG( 2,1 )*BXS(  3, 2 ) -BIG! 3,1 )*BI6( 2,2  > )/DET 
EBIGI 2, 3)  *  -( CIG( 1, 1 )*6IG( 3,2 >-BIG( 3,1 )*8Z6( 1,2) >/D£T 
EBIGI 3,3)  3  ( BIG( 1 , 1 )*B1G( 2 ,2 )-BIG( 2,1 )*BIG( 1,2) )/DET 

FORM  SECOND  DERIVATIVE  MATRIX 

00  131  1=1.2 
DO  131  J=1.3 
EXS(J.I)  =0.00 
00  131  K=1,NEL 

EXS(J.I)  =  EXS(J.I)  ♦  X(I,K)*ESHP(J,K) 

FORM  JACOBIAN  MATRIX  INVERSE 

SXIl.l)  =  XS( 2 ,2 )/XSJ 
SX(  2,2)  =  XS(1,1)/XSJ 
SX( 1 , 2  I  =  -XSU,2)/XSJ 
SX(2,1 1  =  *XS(2,1 l/XSJ 

FORM  GLOBAL  SECOND  DERIVATIVES 


DO  132  1=1. NEL  00000870 

tempi  i  i  =  eskp(i.i)  oooooaao 

TEMPI  2)  3  ESHPI 2,1)  00000890 

tempi  3 )  =  eshpis.d  00000900 

DO  133  J=1 >3  00000910 

ESHPIJ.I)  3  0.00  00000920 

00  134  K=l>3  00000930 

ESHPIJ.I)  =  ESHPIJ.I)  ♦  EBIGI J,K)*(TEMP(K)-  (EXS(K,1>»ISX(1,1>*  Q00C0940 
1  SHP( 1 ,1 )*SX( 1 ,2 )*SHP< 2>I))1~ ( EXSI K,2  )*( SXI 2,1  )*5HP( 1,1 )*SX( 2,2 )*  00000950 


00000590 

00000590 

00000600 

00000610 

00000620 

00000630 

00000640 

00000650 

00000660 

00000670 

00000680 

00300690 

00000700 

00000710 

00000720 

00000730 

00000740 

00000750 

00000760 

00000770 

00000780 

00000790 

00000800 

00000810 

00000820 

00000830 

00000840 

00000850 

00000860 

00000870 

00000380 

00000890 

00000900 

00000910 

00000920 

00000930 


2  SHP(2,I)))) 
continue 
continue 

CONTINUE 

RETURN 

END 

SUBROUTINE  PFORMI  UL  ,  XL  ,  Tt  ,  LD  ,  P  ,  S  ,  IE  , 

1  X  ,  IX  ,  F  ,  T  ,  JDIAG  ,  B  >  A  ,  C  ,NDF, 

2  N0M,NEN1,NST,ISW.U,UD,AFL,BFL.CFL,DFL) 

COMPUTE  ELEMENT  ARRAYS  AND  ASSEMBLE  GLOBAL  ARRAYS 


C*************** ******** 


00000960 

00000970 

00000*80 

00000990 

00001000 

QOOOlOiO 

00000010 

00000020 

00000030 

00000040 

00000050 

*00000060 


101 


ooo  noon  on  uinoo  ooo 


.  i 


£**«#****•»«**«*****«***#  PFORM  #**«#»*»»*»»«*»»#*  00  OCO  70 

>#**#***#*#  <iiH>«a»«*iHH>«»«k,>iMHt«»«ii<HHH>«*«<,»«O0033C30 

C  00000090 

IMPLICIT  REAL*8l A-H.O-Z)  00000100 

LOGICAL  AFL.BFL.CFL.OFL  00000110 

COMMON  /CD AT A/  0  >HEA0l 20 1 tNUMNP ,NUMEL,NUMMAT ,NEN,NEQ,IPR  00000120 

COMMON  /EL0ATA/  DM,N,MA,MCT,IEL,NEL  00000130 

COMMON  /PR LOO/  PROP  00000140 

COMMON  /FVISC/  K2  00000145 

|  COMMON  /TAYLR/  ESIGK 4 ,2 ,50 ) ,ESIG2( 4, 2 ,50 ) ,E5IG3( 4,2,50 ) ,  00000150 

1  TY( 4,2,50 ) ,ELAS1( 4,2,50 ) ,ELAS2( 4,2,50 1  ,ELA33( 4,2,50 ) ,  00000160 

2  BOSIGI 4,2 ,50 ) 

DIMENSION  XL! N0M.1 ) , LDINDF ,1 ) ,P( 1 ) ,SINST, 1 ) , IEI 1 ) ,D( 30.1 ) , ID( NOF 
1 1 , X( NDM, 1 ) , IX( NEN1 ,1 ) , FI NOF , 1 ) , JDIAGI 1 ) >B( 1 ) , A( 1 ) ,C(  1 ) ,UL(N0F ■ 1 1 
2  >TL( 1 ),T( 1 ),U( 1 I,U0(NDF,1) 

C 

IF( I K2. LE .  2 -OR-K2 .  EQ.5  >  .OR .  I ISU. LE.4 ) .OR . (NOF -GE .4 ) )  60  TO  102 
SET  ITERATION  PARAMETERS  FOR  FLUID  VISCOELASTICITY 


NSTEP  =  0 
TOL1  =  l.E+1 

BEGIN  VISCOELASTIC  ITERATION:  LOOP  ON  ELEMENTS 
IEL  *  0 

DO  101  N  *  l.NUMEL 

CALCULATE  ELAS  WITHIN  ELEMENTS  USING  CENTRAL  DIFFERENCES; 
THESE  WILL  BE  USED  FOR  BOUNDARY  ELEMENTS 


GAUSS  POINT  1 

AA  =  YY(4,2,NJ 
E3  =  YY( 2 , 2,N ) 
CC  =  YYI2.1.NJ 
DD  =  YYI4.1.N) 
ELAS1( 1 ,1 ,N )  * 
1  -BOSIGU.l.N) 
ELAS1( 1 ,2 ,N )  = 
1  -BOSIGI 1,1, N> 
ELAS2(1,1,N)  = 
1  -BOSIGU.2.N) 
ELAS2( 1,2 ,N)  * 
1  -BOSIGI 1, 2, N> 
ELAS3( 1,1 ,N)  * 
1  -BOSIGt 1 >  3,N ) 
ELAS3( 1 ,2 ,N)  * 
1  -BOSIGI 1 , 3,N ) 

GAUSS  POINT  4 


-XI 2,1X1 1,N) 
-XI 2 ,IX( 1,N) 
-XI 1,1X1 1 ,NI 
-XI 1,1X1 1 ,N) 
1 1  ESIGK  2,1 
)*BB)/(CC»AA 
( I  ESIG11 4,1 
)*DO)/lCC*AA 
( (ESIG2I 2,1 
)«EB)/ICC*AA 
( I ESIG2I 4 ,1 
)«DO )/( CC*AA' 
1 1 ESIG31 2,1 
)*8B)/(CC»AA 
1 1  ES1G3(4,1 
>*OD  )/(CC*AA 


) 

1 

l 

I 

,N)-BOSIGI 1,1,N) )*AA- 
BB*DD ) 

,N)-BOSIG( 1,1, N) 1«CC- 
BB*D0 ) 

,N ) -BOSIGI 1,2, N) )»AA- 
-BB*DD ) 

,N ) -BOSIGI 1 ,2 ,N) )»CC- 
,-BB*DD ) 

,N)-BOSIG( 1,3,N) )«AA- 
-BB*CO ) 

,N)-BOSIGl 1 >3»N) )*CC- 
-BB*DO ) 


I  ESIGK  4,1  ,N ) 
I E5IG1I 2,1, N) 
(ESIG2(4,1,N> 
I ESIG2I 2,1,N) 
I  ESIGK  4,1  ,N ) 
(ES163I2,1,N) 


AA  :  XI 2, IXI4, N) ) — YYI 1 ,2 ,N ) 

BB  =  YYI 3,2, NI-XI 2,1X1 4,N) ) 

CC  =  YYI 3,1 ,N >— XI 1 ,IX(4,N ) I 
DO  s  XI 1 ■ IXI4, N I )-YY( 1,1,N) 

ELAS1I4 , 1  ,N  I  *  1 1  ESIGK  3,1  ,N)-BOSIG( 4,1  ,N)  )*AA-( BOSIGI 4, 1,N) 
1  -ESIG11 1 , 1 ,N )  >*BB  „  I CC»AA-BB»DD ) 

ELAS114.2.N)  s  1 1  BOSIGI  4,1  ,N ) -ESIGK 1 >1 ,N )  )*CC-(  ESIGK  3,1  ,N) 
1  -BOSIGI 4 , 1 ,N ) )«0D )/l CC»AA-BB»DD ) 


,100000170 

00000180 

00000190 

00000200 

00000210 

00000220 

0C000230 

00000240 

00000250 

00000260 

00000230 

00000290 

00000330 

00000310 

00000320 

00000010 

00000020 

00000030 

00000040 

00000050 

00000060 

00000070 

0C000080 

00000090 

00000100 

00000110 

00000120 

00000130 

00000140 

00000150 

00000160 

00000170 

00000180 

00000190 

00000200 

00000210 

00000220 

00000230 

00000240 

00000250 

00000260 

00000270 

00000280 

00000290 

00000300 

00000310 

00000320 

00000330 


* 

( 


! 


102 


ELAS2(4,1,N)  =  (IESIG2<3,1,N)-B0SIG<4,2,N))*AA-IB0SIG(4,2,N) 
1  -ESIG2I  1 , 1  ,N )  )«E3 )/( CC*AA-B3*DD I 
ELAB:i4,2,N)  a  ( 1 BOSIGI 4, 2 >11  I-ESIG2I  1 >1 ,N I  )*CC-(  ESXG2I  3,1 ,11 1 
1  -EC3IG<  4 , 2 ,N )  )*OD  1/1  CC*AA-B3*DD  t 
ELAS3I 4,1>N)  =  ( I ESIG3I 3.1 ,N)-B0SIG< 4.3.N) )«AA-(B0SIGI4,3.N) 
1  -ESIG3I 1 , 1  ,N )  )*EB  >/( CC*AA-B3*DD ) 

ELAS31 4,2  >N )  a  ( (B0SIG(4,3,N)-ESIG3(1,1,N> >*CC-( ESIG3I 3,1 ,N) 
1  -BOSIGI 4 , 3,N ) 1*00 )/( CC*AA-B3*0D I 

GAUSS  POINT  3 


AA  a  XI 2  > IXI 3>N ) l-YYI 2 • 2 ,N) 

BB  a  XI 2,1X1 3>N) l-YYI 4,2, N) 

CC  a  XII, IXI3, N)I-YY(4,1,N1 
OD  a  XI 1,1X1 3, N) l-YYI 2  >1,N) 

ELASXI  3,1,N)  =  1 I BOSIGI 3,1 ,NI-ESIG1(4,1,N> l*AA 
1  -ESIG1 1 2 , 1 , N  I ) «B3 )/l CC*A A-B3«0D I 
ELASXI 3,2,N)  a  (( BOSIGI 3,1, NI-ESIG1I 2,1, N) )*CC 
1  -ESXGlt  4 , 1 ,N I 1*03  »/( CC«AA-BB*DD ) 

ELAS2I 3,1 >N)  a  ( (BOSIGI 3,2,NI-ESXG2I4,1,N) )«AA 
1  -ESXG2I 2,1,N) )*CB I/I CC*AA-BB*DD I 
EUS2I3.2.NI  a  ( ( BOSIGI  3 , 2,N  I-ESIG21 2 , 1  ,N )  )«CC 
1  -ESXG21 4 , 1  ,N I  I*f00  I/I CCt<AA-B0*.0D  I 
ELAS31 3,1  ,N)  a  ((BOSIGI  3, 3,NI-ESXG3(4,1,NII»AA' 
1  -ESIG3I 2 , 1 ,N 1 )«DB 1/1 CC»AA-B3*D0 I 
ELAS3I 3,2 ,N I  a  ( (BOSIGI 3,3, N1-ESIG3I 2,1, N) )«CC 
1  -ESXG31 4 , 1 ,N 1 1*00 1/1 CC»AA-BB»OD l 


-( BOSIGI 3,1, N) 
-IB0SIG(3,1,NI 
-(BOSIGI 3, 2, Ml 
-(BOSIGI 3, 2, Nl 
-IB0SIG(3.3,NI 
-I BOSIGI 3,3, Nl 


C  GAUSS  POINT  2 
C 

AA  =  YYI 3,2 ,N1-X( 2 , IXI 2 ,N 1 1 
EB  a  X( 2 ,IX( 2,N) l-YYI 1 ,2 ,NI 
CC  =  XI  1,IX(  2,N1  l-YYI  1 ,1  ,N! 

DO  =  YYI  3, 1  ,N l-XI  1 , IXI  2 ,N  1 1 
ELAS1(2,1,NI  a  I (BOSIGI 2,1 ,NI-ESIGI( 1 ,1 ,h 
1  -BOSIGI 2 . 1 ,N I l«B3 I/I CC»AA-BB*DD I 
ELAS1I 2,2,N)  a  ( I ESIG1I 3,1 ,N ) -BOSIGI 2 ,1,K 
1  -ES1G11 1 , 1 ,N 1 1*00  I/I CC*AA-B3*0D  I 
ELAS2(2,1,N)  a  ( (BOSIGI 2,2, N I-ESIG2I 1,1, K 
1  -BOSIGI 2 , 2 ,N I l*BB )/( CC*AA-B8*0D  I 
ELAS2(2,2,N)  a  ( I ESIG2 1 3,1 ,N ) -BOSIGI 2 ,2 
1  -ESIG21 1 ,1 ,H 1 1*00  I/I CC*AA-BB*00 1 
ELAS3I 2 ,1 ,N)  a  ( (BOSIGI 2 ,3,N I-ESIG31 1 ,1>P 
1  -BOSIGI 2 , 3 ,N 1 i*BB I/I CC*A A-BB*00 1 
ELAS3I 2 ,2 ,NI  a  ( I ESIG31 3, 1 ,N J-BOSIGI 2 , 3,h 
1  -ESIG3I 1 ,1 ,N I 1*00 l/l CC*AA-B3*00 1 
C 

C  REPLACE  ELAS  FOR  INTERIOR  ELEMENTS 

C 

DO  91  IDEX  a  l.NUMEL 
00  92  JDEX  =  l.NUMEL 
C 

C  GAUSS  POINT  1 


IFKIXll.Nl.NE. IXI4, IDEXII.OR.  11X12, NI.NE. IXI3, 10EX1IIG0  TO  10 
XF( (IXI l.NI  .NE.IXI2, JDEX) I . OR. I IXI 4, N I. NE. IXI 3, JDEX 1 1 1G0  TO  10 
AA  a  YY(4,2,N)-YY(4,2,I0EX) 

BB  a  YYI 2 ,2,N l-YYI 2,2 , JDEX) 

CC  a  YYI 2 >1,N l-YYI 2,1, JDEX) 

DO  a  YYI 4, 1,N l-YYI 4,1, IDEX I 


A-l ESIG1I 3,1 ,N) 
D-IBOSIGI 2,1,N> 
A-(ESIG2(3,1,N) 
:-IB0SIG(2,2,N) 
A-IESIG3I 3,1,N) 


C-(B0SI6(2,3,N) 


00000340 

00000350 

03C00360 

cocooiro 

000003B0 

00000390 

00000400 

00000410 

00000420 

00000430 

00000440 

00000450 

00000460 

00000470 

00000480 

00000490 

00000500 

00000510 

00000520 

00000530 

00000540 

00000550 

OOOOOESO 

00000570 

00000530 

00000590 

00090600 

00000610 

00000620 

00000630 

00000640 

00000650 

00000660 

000C0670 

00000680 

00000690 

00000700 

00000710 

00000720 

00000730 

00000740 

00000750 

00000760 

00000770 

00000780 

00000790 

00000600 

00000610 

00000820 

00000330 

00000340 

00000350 

00009360 

00000370 

00000380 

00000390 

00000400 

00009410 

00000420 

00000430 


103 


c 

c 

c 

10 


ELASH  1 » 1  iN )  =  ( < ESIGH  2 ,1,N  I -ESIGH  2 , 1 ,  JOEX )  )»AA-(  ESIGH  4,1  ,N  I 
1  -ESIG1<4,1,I0EX>I*S3)/(CC*AA-E3-DD) 

ELASH  1,2,11)  S  (<  ESIGH 4. 1,N)-ESIG1<4,1, IDEX ))#CC-<  ESIGlt  2, 1,N> 
1  -ESIGn2,l.J0EX))*0D)/<CC»AA-ED''0D> 

ELAS21 1 ,1 . N )  s  ( < ESIG2( 2 ,1 ,N  I-ESIG21 2,1, JDEX) )*AA-< ESIG2I 4, 1 ,N ) 
1  -ESIG2I4 , 1 , IBEX ) )*BS )/( CC»AA-E3»DD ) 

ELAS2(1,2,N>  =  <(ESIG2(4,1,N)-ESIG2(4,1,IDEX))«CC-<ESIG2<2,1,N) 
1  -ESIG2C 2 ,1 , JOEX )  1»C0 )/( CC«AA-BB*00 ) 

ELAS3(  1 ,1  >N)  *  (  ( ESIG3I 2 >1  >N)-ES1G3(  2 >1, JDEX) )*AA-(  ESIGH  4,1  ,N) 
1  -ESIG31 4,1.  IOEX I )  <*BB )/( CC«*AA-BB»00  ) 

ELAS3(1,2,N)  =  (( ESIG3< 4,1 ,N)-ESIG3( 4,1, IOEX) )*CC-< ES1G3I 2,1, N) 
1  -ESIG3I 2 , 1 , JDEX >  )«DD  >/< CC«AA-BB*0D ) 

GAUSS  POINT  4 

IF((IX(1,N).NE.XX<2, IOEX)). Oft. (IX(4,N).NE.IX(3,XDEX))>G0  TO  20 


IF(  ( IX( 3,N) .NE. IX( 2, JOEX  ) ) .OR. 
AA  s  YY(1,2,J0EX)-YY(1,2,N> 
2,N)-YY( 3, 2, IOEX) 

1 ,N)-YY( 3,1, IOEX ) 

1 > JOEX )-YY( 1,1, N) 


( IX(4,N).NE. IX(1, JDEX) ) )60  TO  20 


C 

C 

C 

20 


BB  =  YY( 5 
CC  *  YY( 3 
CO  =  YY( 1 
E  LASH  4,1 
1  -ESIGH1 
ELAS1(4,2 
1  -ESIGH3 
ELAS2I4.1 
1  -ESIG2C 1 
ELAS2(4.2 
1  -ESIG2(3 
ELASH  4,1 
1  -ESIG3( 1 
ELAS3(4,2 
1  -ESIG3(3 


00000440 

00000450 

oocoe-iio 

00000470 

00000480 

00000490 

00000500 

00000510 

00000520 

00000530 

00000540 

00000550 

00000560 

00000570 

00000580 

00030590 

00000600 

00000610 

00000620 

00000630 

00030640 


HI  =  ( ( ESIGH  3,1,N)-ESIG1( 3,1 , IOEX) )*AA-( ESIG11 1,1, JDEX 100000650 
1  ,N  I  )#C3 1/1  CC'AA-BB’^OO  I  00000660 
N)  s  ( ( ESIGH  1 ,1 ,  JOEX ) -ESIGH  1 ,1,N)  )»CC-I  ESIGH  3,1  ,N I  000C3670 
1 > IOEX ) I *  DO )/(CC*AA-BB»00 1  00000680 
N1  :  ( ( ESIG2I 3«1,N  >-ESIG2<  3»1»IDEX) )*AA-( ESIG2I 1,1, JDEX >00000690 
1 >N ) )*B8 )/< CC*AA-BB* DO )  00000700 
Nl  3  ( ( ESIG2(  1 ,1 , J0EXI-ESIG21 1,1,N)  )«CC-(  ESXG21 3,1,11)  00000710 
1 , IOEX ) )*00 )/( CC*AA-BB»DO 1  00000720 
Nl  =  ( 1 ESIGH  3,1, N)-ESIG3t 3,1, IDEX Y)*AA-( ESIG31 1,1, JDEX 100000730 


1 ,N ) )«BB )/( CC*AA-EB«00 1 

N)  s  ( ( ESIG3( 1 , 1 , JDEX )-ESIG3( 1 ,1 ,N) )*CC-( ESIG3( 3, 
1 , IOEX 1 1*00 )/<  CC»AA-BS«0D  I 


1»N) 


GAUSS  POINT  3 


00000740 

00000750 

00000760 

00000770 

00000730 

00000790 

OOOOOSOO 

C0000Q10 

00000320 

00000830 

00000340 

00000850 


IF( ( IX( 3,N) .NE.IX( 2 , IOEX 1 1 .OR . ( IX( 4,N) .NE . IX( 1 , IDEX 1 1 1 GO  TO  30 
IF C  € IX( 2,N).NE.IX( 1, JDEX) I .CR.I IXt  2 >N).NE.IX(4, JDEX) I IGO  TO  30 
AA  3  YY( 2 ,2,IDEX)-YY( 2 >2 ,N) 

B3  =  YY( 4,2, JDEX)-YY(4»2 ,N) 

CC  =  YY(4,1,JDEX)-YY(4,1.N) 

DO  s  YY( 2,1, IOEX )-YY( 2,1 ,N) 

E LASH 3,1  >N)  3  (( ESIGH 4,1, JDEX l-ESIGH 4, l,N))»AA-( ESIGH 2,1, IDEXIOOOOOG60 
1  -E5IGH2.1,N))«8B)/<CC»AA-B3*DD>  00000370 

E LASH  3,2  >N)  =  (( ESIGH  2 , 1 , IDEX  l-ESIGH  2 , 1  ,N I  )»CC-(  ESIGH 4, 1 ,  JDEX  I00000CS0 
1  -ESIGH  4, 1,N)  1*00  >/<CC*AA-GB»DD>  00000390 

ELAS2( 3,1 ,N)  3  ( (ESXG2(4, 1, JOEX )-ESIG2(4,l,N) )«AA-t ESIG2( 2,1, IDEX >00000900 
1  -ESIG2( 2 , 1 ,N) l*BB )/( CC*AA-8B*DD I  00000910 

ELAS2( 3,2,N I  3  ( ( ESIG21 2,1, IDEX 1-ESXG2I 2,1,N) )*CC-( E3IG2( 4,1, JDEX IC0000920 
1  -ESIG2t4,l,NI)*00l/(CC*AA-eB*00>  00000930 

ELAS3(3,1,N)  s  I ( ESIG3( 4,1, JDEX )-ESIG3( 4,1 >N) >»AA-(ESIG3( 2 >1,10 EX  100000540 
1  -ESIGH 2 ,1,N) )*BB )/(CC*AA-BB*DQ I  00000950 

ELAS3I  3,2  >N I  s  (( ESIGH  2 ,1 , IOEX  1-ESIG31 2 , 1  ,N )  )*CC-(  ESIG31 4 ,1 ,  JOEX  100000960 


C 

c 

Ic 

30 


1  -ESIG3C 4 , 1  ,N )  )«00  >/( CC«AA-BB«00 I 
GAUSS  POINT  2 

IF t  C IX( 2,N) .NE.IXI 1 , IDEX 1 1 .OR. ( IX( 3,N) .NE .IX(4,IDEX) I IGO  TO  92 
IF( (IX( 1,N).NE.IX(4,JDEX) I .OR. ( IX( 2,N).NE.IX( 3, JOEX) ) IGO  TO  92 
AA  s  YY( 3,2,H)-YY( 3,2, JDEX) 


00000970 

00000930 

00000990 

00001000 

00001010 

00001020 

00001030 


r 

t 

% 

X 

X 

11 

B3  =  YY(1,2,IDEX)-YY(1,2,N) 

00001040 

CC  =  YY(1,1.IDEX)-YY(1,1,N1 

00001050 

DO  =  YY(3,1,N)-YY(3,1,JDEX> 

00001060 

ELASK  2 ,1  ,N )  =  ( (  ESIGK  1 , 1  >  IDEX  l-ESIGK  1,1  >N)  )*AA-(  ESIGK  3>1>N) 

00001070 

1  -ESIGK  3,1.  JDEX)  >»  BB >/< CC*AA-EB*DD  > 

03001080 

ELASK  2 >  2>N)  =  ( ( ESIGK  3, 1>N l-ESIGK  3,1,  JDEX)  )*CC-(  ESIGK  1,1, IDEX  10000 1090  1 

1  -ESIGK  1,1, N)  )*DD )/( CC**AA-BB*DD  ) 

00001100 

ELAS2( 2 ,1 ,N)  =  ( ( ESIG2I 1 ,1 ,IDEX )-ESIG2( 1 ,1,N) >*AA-<  ESIG2I 3,1,N  > 

00001110 

1  -ESIG2( 3, 1 , JDEX ) )*BB )/( CC*AA-BB*00 1 

00001120 

ELAS2I 2,2>N)  =  ( ( ESIG2I 3, 1 ,N )-ESIG2( 3,1, JDEX) >*CC-(ESIG2( 1,1, IDEX >00001130  1 

1  -ESIG2(1,1,NI>»D0)/'(CC«AA-BB*‘DD) 

00001140 

ELAS3( 2,1 ,N)  =  ( ( ESIGK  1,1 , IDEX )-ESIG3t 1,1 »N)  )*AA-(  ESIG3I 3,1,N) 

00001150 

1  -ESIG3( 3,1, JDEX > )*BB )/< CC»AA-68*00 ) 

00001160 

ELAS3( 2,2 ,N)  =  ( ( ESIGK  3,1  ,N)-ESIG3(  3,1 ,  JDEX)  )*CC-( ESIG3I 1,1, IDEX >00001170  i 

1  -ESIG3I 1 , 1 ,N ) >»DD )/< CC»AA-BB*DD ) 

00001180 

92 

CONTINUE 

00001260 

91 

CONTINUE 

00001270 

C 

SET  UP  LOCAL  ARRAYS  FOR  CALCULATING  ESIG(LL,2,N) 

00001280 

DO  5S  1=1, NEN 

00001290 

II  =  IX(I,N) 

00001300 

IF  (II.NE.O)  GO  TO  55 

00001310 

TL(  I )  =  0. 

0C001320 

DO  53  J=1,NDM 

00C01330 

53 

XUJ.I)  =  0. 

00001340 

DO  54  J=1,N0M 

00001350 

UUJ.I)  =  0. 

00001360 

UL( J , I+NEN)  =  0. 

00001370 

54 

LO(J.I)  =  0 

00001330 

GO  TO  58 

00001393 

55 

IIO  =  II*NOF-NOF 

00001400 

NEL  =  I 

00001410 

TL(  I )  =  Till) 

00001420 

DO  56  J=1,NDM 

00001430 

► 

56 

XL(J,I1  =  X(J,II) 

00001440 

DO  57  J=1,NDF 

00001453 

[. 

K  =  IABSt ID( J , II ) ) 

00001460 

■ 

UL(J.I)  =  F(J,II)*PROP 

00001470 

f 

ULtJ.I+NENI  =  UD(  J  ,11 ) 

00001480 

i 

IF  (K.GT.O)  UL(J.I)  =  U(K ) 

00001490 

i 

IF  (DFLJ  K  =  IID  ♦  J 

00001500 

V* 

57 

LO(J.I)  =  K 

00001510 

58 

CONTINUE 

00001520 

C 

FORM  ELEMENT  ARRAY 

00001530 

V 

MA  =  IXtNENl.N) 

00001540 

\ 

IF  (lE(MA).NE.IEL)  MCT  =  0 

00001550 

l 

IEL  =  IE(MA) 

00001560 

CALL  ELMLIB(0(1,MA),UL,XL,IX<1.N>,TL,S,P,NDF,NDM,NST,7) 

00001570 

101 

CONTINUE 

000015S0 

YMAX  =0MAX1(  DABSt  ESIGK  1 , 2 ,1  l-ESIGl (1,1,1))  ,DABS(  ESIG2I  1,2,1) 

00001590 

1  -ESIG2( 1,1,1) ) ,DABS( ESIG3( 1,2,1 >-ESIG3( 1,1,1))) 

00031600 

DO  93  1=1 iNUMEL 

00001610 

f  . 

DO  93  J=1 ,4 

00001620 

XMAX  =0MAX1  ( DABS!  ESIGK  J ,  2 , 1  )-ESIGl ( J ,  1 , 1 ) ) , DABSt  ESIG2 ( J , 2 , 1 ) 

00001630 

1  -ESIG2I J, 1,1)), DABSt ESIG3< J ,2 ,1 >-ESIG3( J , 1,1))) 

00001631 

93 

IF  (XMAX.GT.YMAX)  YMAX=XMAX 

00001632 

F 

IF  ( YMAX. LE. TOLl)  GO  TO  102 

00001633 

r 

NSTEP  =  NSTEP  ♦  1 

DO  90  K=1,NUMEL 

DO  90  J=1 ,4 

ESIGK  J  ,1  ,K )  =  ESIGK  J, 2, K) 

00001634 

► 

II 

i  i 

105 

1-  _ 

90 


..  2 


94 


95 


C 

102 


C 


103 


104 

105 


106 


107 

10S 

C 


C 


110 

1000 

1010 

1020 


C 


ES!G2(J,1.K)  *  ESIG2U.2.K) 

ES1G3U.1.K)  *  ESIG3(J,2,K) 

IF  (NSTEP. GE. 10)  GO  TO  102 
IF  ( NSTEP .  EQ . 1 . OR .  NSTEP . EQ . 3 . OR .NSTEP . EQ . 5  - 
1  OB. NSTEP. EQ.9)  GO  TO  94 
GO  TO  5 

WRITE  (6,10001  0, HEAD, TIME, NSTEP 
WRITE  (6,1010) 

00  95  m.NUMEL 

WRITE  (6,1020)  I,((YY( J,1,I),YY( J  ,2 , 1  > ,  ESIGKJ,  2 ,1 ) 

1  ,ESIG2( J,2,I),ESIG3( J,2,I) ),J=1,4) 

GO  TO  5 

LOOP  ON  ELEMENTS:  ELASTIC  ITERATION  COMPLETE 

CONTINUE 

IEL  =  0 

DO  110  N  =  l.NUMEL 
SET  UP  LOCAL  ARRAYS 
00  108  I  =  l.NEN 
II  =  IXC I, N) 

IF  (II.NE.O)  GO  TO  105 
TL(  I )  =  0. 

DO  103  J=1,N0M 
XL(J,I)  s  0. 

DO  104  J  =  1,NDF 
UL(J.I)  =  0. 

UL(J,I+NEN)  =  0. 

LOU,!)  =  0 

GO  TO  108 

IID  =  II«NDF  -  NDF 

NEL  =  I 

TL(  I )  =  T(II) 

00  106  Jsl.NDM 
XL(  J ,  I )  =  X(J, II) 

DO  107  JU.NOF 
K  =  IASS( ID ( J , II ) ) 

UL(J,I)  s  F( J > II )*PROP 
UL( J.I+NEN)  =  UD( J ,11 ) 

IF  (K.GT.O)  UL(J.I)  =  U(K ) 

IF  (DFL )  K  =  IID  ♦  J 
LD(  J ,  I )  =  K 
CONTINUE 

FORM  ELEMENT  ARRAY 
MA  s  IX(NENl.N) 

IF( IE(MA) .NE. IEL)  MCT  =  0 
IEL  =  IE(MA) 

CALL  ELMLIB(0(1,MA),UL,XL,IX(1,N),TL,S,P,N0F,N0M,NST,ISM) 

ADO  TO  TOTAL  ARRAY 

IF( AFL.OR.BFL.OR.CFL)  CALL  AOOSTF( A,B,C,S,P, JDIAG, LO,NST,NEL*NDF, 
1  AFL.8FL.CFL) 

CONTINUE 

FORMAT(A1,20A4,//5X, 'ELASTIC  FLUID  STRESSES  AT  GAUSS  POINTS', 

1  5X, 'TIME', G13. 5, //IX, 'VISCOELASTIC  ITERATION  NUMBER :', 14// ) 
FORMAT) IX, 'ELMT' ,11X, '1-COORO ■ ,6X, ' 2-COORD ' ,20X, 

1  ' ETAU-XX ' , 7X , • ETAU- YY ' , 7X , ' ETAU-XY • // ) 

FCR(1AT(  I5>/4(  10X,2G13.4,13X,3G13.4/)//) 

RETURN 

END 

SUBROUTINE  CMATRX(C,J,SIG,SHP) 


«»«»4*',»#***»*»»¥»*#»**»»»*** 


00001635 

00001636 

00001637 

00001638 

00001639 

00001640 

00001641 

00001642 

00001643 

00001644 

00001649 

00001650 

00001660 

00001670 

00001680 

00001690 

00001700 

00001710 

00001720 

00001730 

00001740 

00001750 

00001760 

00001770 

00001700 

00001790 

00001600 

00001810 

00001811 

00001812 

00001813 

00001814 

00001815 

00001616 

00001317 

00001318 

00001819 

00001820 

00001821 

00001822 

00001823 

00001824 

00001625 

00001826 

00001827 

00001323 

00001829 

00C01830 

00001831 

00001832 

00001333 

00001834 

00001635 

00001849 

00001850 

00000010 

00000020 

00000030 
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C**»******»»*»*r««»»*»*** 

c 

IMPLICIT  REAL»8(A-H,0-Z) 

DIMENSION  C(  1  > >SIG( 7 ) ,SHP( 3,9 ) 

C 

CALL  PZEROIC.6) 

C 

C  ONLY  20  FLOW  TREATED  HERE 
C 

Ctl)  :  2.D0*(SIG(1)*SHP( 1, J )+SIG( 3)*SHP( 2, J ) ) 

C(2)  s  0.00 

C( 3)  3  SIG(2)«SHP(2,J)+SIG(3)*SHP(1,J) 

C(4>  3  0.00 

C(5)  3  2 .00*1 SIG( 2  )*SHP( 2 > J ) +SIG( 3 )»SHPt 1 » J ) ) 

C(6)  3  SIG( 1 )*SHP( 1 , J+SIG( 3 )*SHP( 2 . J ) ) 

RETURN 

ENO 

SU3R0UTINE  FPSIG  (XX.ESIGl , ESIG2 .ESIG3.SIG.ITYPE ,NOF ) 

C 

C************************  FPSIG  «***»««»»»»»**»***•»*«»«*** 
C»*****«*«**MM*«M**«*«*»*  »**#*»*»*«*«**»****#«**»*** 

c 

IMPLICIT  REAL«8(A-H,0-Z) 

DIMENSION  XX(1),SIG(1> 

COMMON  /CDATA/  0,HEA0( 20 ) >NUMNP,NUMEL,NUMMAT,NEN,NEQ,IPR 
COMMON  /ELOATA/  OM,N,MA,MOT,IEL,NEL 
COMMON  /TDATA/  TIME,0T,C1,C2,C3,C4,C5 
COMMON  /FVISC/  K2 
C 

GO  TO  ($1,52,53,54),  ITYPE 
C 

C  PLANE  FLOW 

C 

51  MOT=MOT-l 

IF  (K2.LE.2)  GO  TO  509 
A  3  SIG(1)  ♦  ESIG1 
B  =  SIG( 2 )  ♦  ESIG2 
C  =  SIG( 3)  ♦  ESIG3 
509  IF  (MOT.GT.O)  GO  TO  510 

IF  (NOF.LT.4)  WRITE  (6,5000)  0, HEAD, TIME 
5000  FORMAT  ( A1 , 20 A4 , //5X , ' F LUIO  VISCOUS  STRESSES  AT  GAUSS  POINTS:1, 

1  5X, "TIME ' ,G13.5, 

2  //IX.'ELMT  MATL' ,6X, 'l*COORO' ,6X, '2-COORO' ,5X, 

3  'FRESSURE' ,7X, ’TAU-XX* ,7X, 'TAU-YY' ,7X, 'TAU-XY'/) 

IF  (NDF . LT.4.AND.K2.GE.3)  WRITE  (6,5010) 

5010  FORMAT  (//5X,‘TOTAL  VISCOUS  AND  ELASTIC  STRESSES  ATGAUSS  POINTS: 
I1/) 

IF  (N0F.GE.4)  WRITE  (6,5020)  0, HEAD, TIME 
5020  FORMAT  ( A1 .20A4.//5X, 'TOTAL  VISCOUS  AND  ELASTIC  STRESSES  AT 

1  GAUSS  POINTS: 1 ,5X> 'TIME* ,G13.5> 

2  //IX.'ELMT  MATL' >6X> '1-COORD' >6X» '2-COORD' »5X> 

3  'PRESSURE' ,7X, "TAU-XX" ,7X, 'TAU-YY* ,7X, 'TAU-XY'/) 

IF  (NDF.GE.4)  GO  TO  509 

IF  (K2.GE.3)  M0T=19 
IF  (K2.GE.3)  GO  TO  508 
,  MOT  s  50 

508  CONTINUE 


00000040 

00000050 

00000060 

00000070 

00000080 

00000090 

00000100 

00000110 

00000120 

00000130 

00000140 

00000150 

00000160 

00000170 

00000180 

00000190 

00000200 

00000210 

00000010 

00000020 

00000030 

00000040 

00000050 

00000060 

00000070 

00000080 

00000090 

00000100 

00000110 

00000120 

00000130 

00000140 

00000150 

00000160 

00000170 

00000240 

00000250 

00000260 

00000270 

00000180 

00000190 

00000200 

00000210 

00000220 

00000230 

00000280 

00000290 

0C000300 

00000301 

00000302 

00000305 

00000306 

00000307 

00000308 

00000310 

00000320 

00000330 

00000350 


510  IF  (NDF.LT.4)  WRITE  (6, 5001)N,MA,XX(l),XX(2),SIG(7),(SIG(I), 1=1, 3)00000360 


IF  (NDF.GE.4)  WRITE  (6,5009)N,MA,XX( 1 ) ,XX( 2 ) 


00000363 


.  1 


00000366 

00000370 

00000300 

00000390 

00000400 

00000410 

00000420 

00000430 

00000440 

00000450 

00000460 

00000470 

00000430 

00000490 

00000500 

00000510 

00000520 

00000530 

00000540 

03000550 

00000560 

00000570 

00000560 

00000590 

00000600 


l\ 


5009  FORMAT  (2X5,2613.4) 

IF  (K2.GE.3)  WRITE  (6,5011)  SIG(  5  )  ,SIG( 6  ) , A,B,C 

5001  FORMAT  (2X5.6G13.4) 

5011  FORMAT  (IX, 'RE  =  ‘,G13.4.'WS  =  ' ,G13.4,13X,3G13.4) 

RETURN 

C 

52  RETURN 
C 

C  AXISYMMETRIC  FLOW 

C 

53  M0T=I10T-1 

IF  (MOT.GT.O)  GO  TO  530 
WRITE  (6,5002)  O.HEAO.T1ME 

5002  FORMAT  ( A1,20A4,//5X, ' FLUID  STRESSES  AT  GAUSS  POINTS:', 

1  5X, 'TIME* ,G13.5, 

2  //lX.'ELMT  MATL',6X, '1-COORD*, 6X, '2-COORO* ,5X, 

3  ' PRESSURE ' , 7X , ' TAU-RR ' . 7X , ' T AU-ZZ  * , 7X , ' TAU-TT • , 7X , ' TAU-RZ • / ) 
MOT  =  50 

C 

530  WRITE  (6.5003)  N,MA,XX( 1),XX( 2),SIG( 7),(SIG( I), 1=1, 4) 

5003  FORMAT  (2I5.7G13.4) 

RETURN 

C 

54  RETURN 
END 
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APPENDIX  4 


Input  Data  Set  Listings 

1.  Run  1  -  Linear  Cross  Channel  Flow 
18-9  Node  Elements 

2.  Run  3  -  Linear  Cross  Channel  Flow 
72-8  Node  Elements 

3.  Run  4  -  Convection  (Re  =  0.4)  Cross  Channel  Flow 
18-9  Node  Elements 

4.  Run  6  -  Viscoelastic  (Ws  =  0.02)  Cross  Channel  Flow 
18-9  Node  Elements 

5.  Run  13  -  Viscoelastic  (Ws  =  0.001)  Entry  Flow 
24-9  Node  Elements 

6.  Run  20  -  Linear  Entry  Flow  Fully  Developed  Boundary 
Conditions  24-9  Node  Elements 

(Note:  Run  20  is  not  listed  in  Table  1) 
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BRC4066  (FOREGROUND):  OUTFUT  FROM  TSO  X PRINT 
AT  13:42:06  ON  12/05/80  -  ERC4066.TEST.DAT A 


Input  Dataset  Run  No.  1 


i 


FEAP 

91 

COOR 


CROSS-CHANNEL  FLOW  -  NEWTONIAN  (TEST  7) 
18  1  2  2  9  0 


1 

7 

000 

0D0 

85 

0 

200 

000 

2 

7 

000.166666700 

66 

0 

200.166666700 

3 

7 

000.333333300 

87 

0 

200.333333300 

4 

7 

000 

.500 

88 

0 

200 

.500 

5 

7 

000.666666700 

89 

0 

2D0. 666666 700 

6 

7 

000.833333300 

90 

0 

200.833333300 

7 

7 

000 

100 

91 

0 

200 

100 

ELEM 


1 

1 

1 

15 

17 

3 

8 

16 

10 

2 

9 

14 

7 

1 

3 

17 

19 

5 

10 

18 

12 

4 

11 

14 

13 

1 

5 

19 

21 

7 

12 

20 

14 

6 

13 

14 

MATE 

1  5  NINE-NODE  LAGRAGIAN  PENALTY  ELEMENT 

10  2  1  1.0 

2  .1000+009  .7900+003  .0000 


BOUN 

1  7-1-1 

85  0  1  1 

2  1-1-1 

6  0  11 

86  1  -1  -1 

91  0  1  1 

7  7-1-1 

89  0  1  1 

FORC 

7  7  -102  000 

91  0  -102  000 


ENT) 

_MACR 
UTAN 
FORM 
SOLV 
OISP 
STRE 
RE  AC 
END 

_ STOP 


00000010 

00000020 

00000030 

00000040 

00000050 

00000060 

00000070 

ooooooso 

00000090 

00000100 

00000110 

00000120 

00000130 

00000140 

00000150 

00000160 

00000170 

OOOOOISO 

000C0190 

00000200 

00000210 

00000220 

00000230 

00000240 

00000250 

00000260 

00000270 

00000275 

00000200 

00000290 

00000500 

00000310 

00000320 

00000330 

00000340 

00000350 

00000360 

00000370 

00D003S0 

00000390 

00000400 

00000430 

00000490 

00000500 

00000510 

00000520 

00000530 

00000540 

00000550 

00000560 

00000570 

00000580 


I 
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BRC4066  ( FOREGROUND  ) :  OUTFUT  FROM  TSO  XPRINT 


Input  Dataset  Run  No.  3 


AT  13:42:48  ON  12/05/80  -  BRC4066 . TEST2 . DATA 


FEAP  CROSS-CHANNEL  FLOW— LINEAR  NEWT0NIAN//72  ELEMENTS 


00000010 


ELEM 


21  23 


14 


22  15 


20 


00000540 

00000550 

00000560 


253 

DOT 

72 

1  2 

2  8 

0 

00000020 

00000030  i 

1 

1 

000 

OOO 

00000040  > 

13 

0 

000 

100 

00000050  i 

14 

1. 

083333300 

000 

00000060  ! 

20 

0. 

083333300 

100 

00000070  j 

21 

1. 

166666700 

000 

ooooooso  f 

33 

0. 

166666700 

100 

00000090  < 

34 

1 

.2500 

000 

00000100  1 

40 

0 

.2500 

100 

00000110  1 

41 

1. 

333333300 

000 

00000120  i 

53 

0. 

333333300 

100 

00000130  ! 

54 

1. 

416666 7D0 

000 

00000140 

60 

0. 

416666700 

100 

00000150 

61 

1 

.SDO 

0D0 

00000160 

73 

0 

.500 

100 

00000170 

74 

1. 

533333300 

000 

00000180 

80 

0. 

533333300 

100 

00000190 

81 

1. 

666666700 

000 

00000200 

93 

0. 

666666700 

IDO 

00000210 

94 

1 

.7500 

000 

00000220 

100 

0 

.7500 

100 

00000230 

101 

1. 

833333300 

000 

00000240 

113 

0. 

833333300 

100 

00000250 

114 

1. 

916666700 

000 

00000260 

120 

0. 

916666700 

100 

00000270 

121 

1 

100 

000 

00000280 

133 

0 

100 

100 

00000290 

134 

11 

.08333300 

000 

00000300 

140 

01 

.03333300 

100 

0C000310 

141 

11 

.16666700 

coo 

00000320 

153 

01 

.16656700 

IDO 

00000330  j 

154 

1 

1.2500 

OOO 

00000340 

160 

0 

1.2500 

100 

00000350 

161 

11 

.33333300 

OOO 

00000360  ! 

173 

01 

.33333300 

100 

00000370 

174 

11 

.41666700 

000 

00000380 

180 

01 

.41664700 

100 

00000390 

181 

1 

1.500 

000 

00000400 

193 

0 

1.500 

100 

00000410 

194 

11 

.53333300 

000 

00000420 

200 

01 

.53333300 

100 

00000430 

201 

11 

.66666700 

0D0 

00000440 

213 

01 

.6666670 0 

ICO 

00000450 

214 

1 

1 . 7500 

coo 

00000460 

220 

0 

1,7500 

100 

00000470 

221 

11 

.83333300 

ODO 

00000430 

233 

01 

.83333300 

100 

00000490 

234 

11 

.91666700 

ODO 

00000500 

240 

01 

.91666700 

100 

00000510 

241 

1 

200 

ODO 

00000520 

253 

0 

200 

100 

00000530 

111 


i 


13 

1 

3 

23  25 

5  15 

24 

16 

4 

20 

00000570 

25 

1 

5 

25  27 

7  16 

26 

17 

6 

20 

00000580 

37 

1 

7 

27  29 

9  17 

28 

16 

8 

20 

00000590 

49 

1 

9 

29  31 

11  IS 

30 

19 

10 

20 

00000600 

61 

1 

11 

31  33 

13  19 

32 

20 

12 

20 

00000610 

00000620 

HATE 

00000630 

1 

5 

EIGHT-NOTE  SERENDIPITY 

PENALTY  ELEHENT 

00000640 

1 

0 

1 

1 

1.0 

00000650 

2 

.1000*009 

.7900*003 

.0000 

00000660 

00000670 

BOUN 

00000680 

1 

1 

-1 

-1 

00000690 

13 

0 

1 

1 

00000700 

14 

20 

-1 

-1 

00000710 

234 

0 

1 

1 

00000720 

21 

20 

-1 

-1 

00000730 

221 

0 

1 

1 

00000740 

241 

1 

-1 

-1 

00000750 

253 

0 

1 

1 

00000760 

20 

20 

-1 

-1 

00C00770 

240 

0 

1 

1 

00000780 

33 

20 

-1 

-1 

00000790 

233 

0 

1 

1 

OOOOCSOO 

00000810 

FC3C 

00000820 

20 

70 

-1D2 

000 

00000330 

240 

0 

-102 

000 

00000340 

13 

20 

-1D2 

000 

00000350 

253 

0 

-102 

000 

00000360 

00000370 

END 

00000830 

MACR 

00000390 

TANG 

00000900 

FC.7H 

00000910 

SOLV 

00000920 

OISP 

00000930 

STRE 

00000940 

END 

00000950 

STOP 

00000960 

u 
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Input  Dataset  Run  No. 


B3C4C66  (FOREGROUND):  OUTPUT  FROM  TJO  XPRINT 
AT  16:22:13  ON  12/05/80  -  BRC4066 . TEST . DATA 


FEAP  CROSS-CHANNEL  FLOW  -  NEWTONIAN  WITH  CONVECTION 
91  18  1  2  2  9  0 

coca 

1  7  ooo  ooo 

85  0  200  000 

2  7  000.166666700 

86  0  200.166666700 

3  7  000.333333300 

87  0  200.333333300 

4  7  OOO  .500 

88  0  200  .500 

5  7  000.666666700 

89  0  200.666666700 

6  7  000.833333300 

«0  0  200.833333300 

7  7  0D0  100 

91  0  200  100 


ELEM 

1  1 
7  1 

13  l 


1  5  NINE-NODE  LAGRAGIAN  PENALTY  ELEMENT 

10  11  1.0 

2  .1000*009  .7900*003  1.6000 


17 

3 

8 

16 

10 

2 

9 

14 

19 

5 

10 

18 

12 

4 

11 

14 

21 

7 

12 

20 

14 

6 

13 

14 

-1  -1 

1  l 

-1  -1 

1  1 

-1  -1 

1  1 

-1  -1 

1  1 


00000010 

00000020 

00000030 

00000040 

00000050 

00000060 

00000070 

00000080 

00000090 

00000100 

00000110 

00000120 

00000130 

00000140 

00000150 

OCOOOlbO 

00000170 

00000180 

00000190 

0-000200 

00000210 

00000220 

00000230 

00000240 

03000250 

00000260 

00000270 

00000275 

00000200 

00000290 

00000300 

00000310 

00000320 

00000330 

00000340 

00000350 

00000360 

00000370 

00000330 

00000390 

00300400 

00000480 

00003490 

00000500 

00000505 

00000520 

00030530 

00000540 

03000550 

00000560 

00000565 

00000567 

00000570 

00000580 

00000590 

00000595 


STOP 


00000600 

00000610 
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BRC4066  ( FOREGROUND  ) !  OUTPUT  FROM  TSO  XPRINT 

AT  13:50:54  ON  12/07/S0  -  ERC4066. TEST. DATA 


Input  Dataset  Run  No.  6 


FEAP  SQUARE  CAVITY-  OLOROYO  VISCOELASTIC  <RH0=0.  P4=i,  WS= . 01 1 


91 

18 

1 

2  2 

9 

0 

COCR 

1 

7 

000 

000 

85 

0 

200 

000 

2 

7 

000.166666700 

86 

0 

200.166666700 

3 

7 

0D0.3333333D0 

87 

0 

200.333333300 

4 

7 

000 

.500 

88 

0 

200 

.500 

5 

7 

000.666666700 

89 

0 

200.666666700 

6 

7 

000.833333300 

90 

0 

200.833333300 

7 

7 

000 

100 

91 

0 

200 

100 

ELEM 

1 

1 

1 

15  17 

3 

8 

16 

10 

2 

9 

7 

1 

3 

17  19 

5 

10 

18 

12 

4 

11 

13 

1 

5 

19  21 

7 

12 

20 

14 

6 

13 

MATE 

1 

5 

NINE-NODE  LAGRAGIAN 

PENALTY 

ELEMENT 

1 

0 

3 

1 

1.0 

2  .1000*009  .7900*003  .3950*007 

boun 

1  7-1-1 

85  0  1  1 

2  1-1-1 

6  0  11 

86  1  -1  -1 

91  0  1  1 

7  7-1-1 

89  0  1  1 

FORC 

7  7  -1D2  000 

91  0  -102  0D0 


END 

MACR 

“dt 

LOOP 

UTAH 

FORM 

SOLV 

OXSP 

STRE 

TIME 

“NEXT 

OISP 

STRE 

REAC 


1. 

20 


00000010 

00000020 

00000030 

00000040 

00000050 

00000060 

00000070 

00000080 

00000090 

00000100 

00000110 

00000120 

00000130 

00000140 

00000150 

00000160 

00000170 

00000180 

00000190 

00000200 

00000210 

00000220 

00000230 

00000240 

00000250 

00000260 

00000270 

00000275 

00000280 

00000290 

00000300 

00000310 

00000320 

00000330 

00000340 

00000350 

00000360 

00000370 

00000380 

00000390 

00000400 

00000485 

00000490 

00000500 

00000510 

00000520 

00000530 

00000540 

00000550 

00000558 

00000566 

00000575 

00000585 

00000588 

00000591 

00000595 


! 
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00000600 

00000610 


6RC4066  (FOREGROUND):  OUTPUT  FROM  TSO  X PRINT 


Input  Dataset  Run  No.  13 


AT  15:38:57  ON  12/07/60  -  ERC4Q66 . TEST 3 . DATA 


FEAP  ENTRY  FLOW  -  OIDROYD  VISCOELASTIC  (RHO*0.  PA*:, 


121 

24 

1  2 

2  9 

COOR 

1 

1 

000 

000 

9 

0 

000 

100 

10 

1 

.12500 

000 

16 

0 

.12500 

100 

19 

1 

.2500 

ODO 

27 

0 

.2500 

100 

26 

1 

.37500 

000 

36 

0 

.37500 

100 

37 

1 

.500 

000 

45 

0 

.500 

IDO 

46 

1 

.62500 

000 

54 

0 

.62500 

100 

55 

1 

.7500 

ODO 

63 

0 

.7500 

100 

64 

1 

.87500 

ODO 

72 

0 

.87500 

100 

73 

1 

100 

000 

61 

0 

100 

100 

82 

5 

1.12500 

.2500 

117 

0 

200 

.2500 

63 

5 

1.12500 

.37500 

118 

0 

200 

.37500 

84 

5 

1.12500 

.500- 

119 

0 

200 

.500 

85 

5 

1.12500 

.92500 

120 

0 

200 

.92500 

86 

5 

1.12500 

.7500 

121 

0 

200 

.7500 

MS*0.0005> 


ELEM 


1 

1 

1 

19 

21 

3 

10 

20 

12 

2 

11 

16 

5 

1 

3 

21 

23 

5 

12 

22 

14 

4 

13 

16 

9 

1 

5 

23 

25 

7 

14 

24 

19 

9 

15 

16 

13 

1 

7 

25 

27 

9 

16 

29 

18 

6 

17 

16 

17 

l 

75 

87 

89 

77 

62 

86 

84 

76 

83 

0 

16 

1 

87 

97 

99 

69 

92 

98 

94 

88 

93 

10 

21 

1 

77 

69 

91 

79 

64 

90 

66 

76 

65 

0 

22 

1 

69 

99 

101 

91 

94 

100 

99 

90 

95 

10 

MATE 

1  5  NINE-NODE  LAGRANGE  PENALTY  ELEMENT 

10  3  1  1.0 

2  .1000*009  .7900*003  .0000  .79SO*OOS 


BOUN 

1  9-1-1 

73  0  1  1 

2  l  -1  -l 

,  8  0  11 
.  9  9-1-1 

61  0  1  1 

79  0  1  1 

75  0  1  1 


l 


00000020 

00000038 

00000090 

00000050 

00000090 

00000070 

00000880 

00000090 

00000100 

00000110 

00000120 

00000130 

00000190 

00000150 

00000160 

00000170 

00000180 

00000190 

00000200 

00000210 

00000220 

00000230 

00000290 

00000250 

00000290 

00000270 

00000260 

00000290 

00000300 

00000310 

00000320 

00000330 

00000340 

00000350 

00000390 

00000370 

00000380 

00000390 

00000900 

00000410 

00000920 

00000930 

00000440 

00000450 

00000460 

00000470 

00000480 

00000490 

00000500 

00000510 

00000520 

00000530 

00000540 

00000550 

00000590 
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Input  Dataset  Run  No.  20 


BRC4066  ( FOREGROUND ) :  OUTPUT  FROM  TSO  XFRINT 

AT  18=59:50  ON  12/18/80  -  BRC4066.TEST3.DAT A 


FEAP  ENTRY  FLOW  -  OLOROYD  VISCOELASTIC  (RHO=0,  P4=l,  WS*0.0005> 
121  24  l  2  2  9  0 

COOR 


1 

1 

000 

000 

9 

0 

000 

100 

10 

1 

.12500 

000 

16 

0 

.12500 

100 

19 

1 

.2500 

000 

27 

0 

.2500 

100 

26 

1 

.37500 

000 

36 

0 

.37500 

100 

37 

1 

.500 

000 

45 

0 

.500 

100 

46 

1 

.62500 

000 

54 

0 

.62500 

100 

55 

1 

.7500 

000 

63 

0 

.7500 

100 

64 

1 

.87500 

000 

72 

0 

.87500 

100 

73 

1 

100 

000 

81 

0 

100 

100 

82 

5 

1.12500 

.2500 

117 

0 

200 

.2500 

83 

5 

1.12500 

.37500 

118 

0 

200 

.37500 

84 

5 

1.12500 

.500 

119 

0 

200 

.500 

85 

5 

1.12500 

.62500 

120 

0 

200 

.62500 

66 

5 

1.12500 

.7500 

121 

0 

200 

.7500 

ELEN 


1 

1 

1 

19 

21 

3 

5 

1 

3 

21 

23 

5 

9 

1 

5 

23 

25 

7 

13 

1 

7 

25 

27 

9 

17 

1 

75 

87 

89 

77 

16 

1 

67 

97 

99 

89 

21 

1 

77 

89 

91 

79 

22 

1 

89 

99 

101 

91 

MATE 

1 

5 

NINE-NODE  LAGRANGE 

1 

0 

1 

1 

1.0 

2  .1000+009  .7900+003 


10 

20 

12 

2 

11 

16 

12 

22 

14 

4 

13 

18 

14 

24 

16 

6 

15 

18 

16 

26 

IS 

8 

17 

18 

82 

88 

64 

76 

83 

0 

92 

98 

94 

88 

93 

10 

84 

90 

86 

78 

85 

0 

94 

100 

96 

90 

95 

10 

PENALTY  ELEMENT 
00  .7950+008 


BOUN 

1  9-1-1 

73  0  1  1 

75  0  1  1 

79  0  1  1 

2  10-1 

8  0  0  1 

9  9-1-1 

«1  0  1  1 


00000010 

00000020 

00000030 

00000040 

00000050 

00000060 

00000070 

00000060 

00000090 

00000100 

00000110 

00000120 

00000130 

00000140 

00000150 

00000160 

00000170 

00000180 

00000190 

00000200 

00000210 

00000220 

00000230 

00000240 

00000250 

00000260 

00000270 

00000230 

00000290 

00000300 

00000310 

00000320 

00000330 

00000340 

00000350 

00000360 

00000370 

00000380 

00000390 

00000400 

00000410 

00000420 

00000430 

00000440 

00000450 

00000460 

00000470 

00000480 

00000490 

00000500 

00000503 

00000506 

00000510 

00000520 

00000530 

00000540 


74 

80 

82 

1X7 

66 

121 

118 

119 

120 


FORC 


END 

MACS 

UTAN 

FORM 

SOLV 

OISP 

STRE 

REAC 

END 

STOP 


0 

0 

5 

0 

5 

0 

0 

0 

0 


1 

1 

-1 

1 

-1 

1 

0 

0 

0 


1 

1 

-1 

1 

-1 

1 

1 

1 

1 


2  0 

3  0 

4  0 

5  0 

6  0 

7  0 

8  0 


4.202 

2.102 

4.202 

2.102 

4.202 

2.102 

4.202 


000 


000 


00000550 

oeocosco 

00000590 

00000600 

00000610 

00000620 

00000630 

00000640 

00000650 

00000660 

00000670 

00000660 

00000681 

00000682 

00000663 

00000664 

00000685 

00000690 

00000700 

00000710 

00000720 

00000750 

00000760 

00000770 

00000780 

00000790 

00000840 

00000850 

00000660 


APPENDIX  5 


Brief  Review  of  Gyroscope  Theory 

This  Appendix  is  presented  for  the  benefit  of  the  materials 
engineer  who  may  not  be  familiar  with  the  theory  of  gyroscopic 
behavior.  The  discussion  is  taken  entirely  from  Wrigley  et.  al. 
[36].  Figure  18  shows  a  cutaway  of  the  single  degree  of  freedom 
gyroscope  used  in  this  study.  The  normal  assumptions  for  the 
description  of  a  gyro  element  performance  are: 


r\ 


1.  The  rotor  spins  about  an  axis  of  symmetry. 

2.  The  rotor  spins  at  constant  speed. 

3.  Spin  angular  momentum  is  much  greater  than  non-spin 
angular  momentum. 

4.  Center  of  mass  of  the  rotor  and  gyro  element  coincide, 
and  5.  The  rotor  bearing  structure  is  rigid. 

For  a  platform  stabilized  single  degree  of  freedom  gyro, 
these  assumptions  lead  to  the  performance  equation: 


—  +  °g  i  +  V  =  Hsf"lA  '  “=md  -  9“SRA 


3  dt* 


da), 


-  I 


OA 


g  dt 


For  integrating  gyros,  a  restraining  torsional  spring  is 
eliminated,  hence  kg  =  0  and  the  performance  equation  becomes: 


or 


d20  .  d0  _  Hs  f 
g  dt7  dt  Cg  [' 

3t(Tg  H  +  9)  ~  3t( 


WIA  "  “cmd  “  0WSRA 


d“0A 
g  dt 


IA  “  “cmd  ’  0WSRA 
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•  ^ «)  • 


T 

g  OA 


TT  ~~Z 


Therefore: 


T 

g 


e 


Where 


0), 


OA 


0  =  Output  axis  rotation 
oj^  5  Input  axis  angular  rate 

“crad  H  Commanded  output  axis  angular  rate 
“sra  5  sPin  reference  axis  angular  rate 
H  =  Rotor  angular  momentum 

S 

t  =  I  /c  =  time  constant 
y  y  y 

I  =  gyro  output  axis  effective  moment  of  inertia 
Cg  =  float  damping  coefficient 
U^Mqa  j  =  Uncertain  torgue  about  output  axis 


Assuming  0  &  Tg  <<1  &  u>IA  =  “cmd'  eguation  becomes 


Tg  f  +  8  -  ^/"("oa)  at 

This  equation  shows  that  the  gyro  drift  uncertainty  is  a  first 
order  response  to  the  time  integral  of  the  uncertainty  torques 
about  the  output  axis. 


Alternately  expressing  the  equation  in  terms  of  drift  rate. 


r  d“0A  , 

Tg  ~3t”  +  woa 


u(moa) 


Hence,  any  source  of  uncertain  torque  of  the  torque  summing 
member  about  the  output  axis  is  a  contributor  to  the  possible 
inaccuracy  of  the  gyro  element. 
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The  most  common  sources  of  these  uncertainties  are 
gimbal  friction  and  mass  unbalance. 

These  are  factors  very  sensitive  to  the  material  state 
and  processing  variables.  It  is  for  this  reaso..  that  a 
rational  method  of  selecting  injection  molding  parameters  is 
required.  A  brief  example  of  this  is  presented. 

Prior  to  introduction  into  service,  the  molded  gyro  is 
balanced.  Remaining  unbalance  can  be  nullified  by  compensa¬ 
tion  in  the  feedback  loop  of  the  control  system.  However, 
from  the  drift  rate  equation  we  can  see  that  for  a  step 
acceleration  the  steady  state  (t  -*■  °°)  drift  rate,  due  to 
torque  uncertainties  caused  by  variations  to  the  balance,  is 


U). 


OA 


-  PVeagTg 


s.s. 


Where  p  is  the  mass  density,  V  is  the  effective  volume 
of  unbalance,  e  is  the  amount  of  mass  eccentricity,  and  a  is 
the  step  acceleration  in  g's. 

Taking  typical  values: 

c  =  20588  <&2£zsec 

g  cm 

t  =  0.0017  sec 

g 


“oaI  =  l0//hr  =  4.8  x  10~®  rad/sec 
's.s. 

p  =  1.6  gm/cm^  (Polypheneline  Sulfide) 
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We  obtain: 


_  0.0363  _ 4 

v  e  *  - - —  cm  . 


For  an  acceleration  of  lOgs  then  we  get: 


Ve  »  0.00363  cm4 

which  defines  the  bounds  of  mass  imbalance  which  can  be  tolerated, 
for  the  specified  performance,  due  to  long  term  materials  be¬ 
havior  (creep  relaxation,  non-uniform  thermal  strain,  etc.). 
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